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A  Study  of  the  Free  Energy  of  the  Lenz-Ising  Model 
Using  the  Cluster  Variation  Method  of  Morita 

CHAPTER  I 

CLUSTER  VARIATION  METHOD 
Objectives  of  This  Work 
The  purpose  of  this  work  can  be  described  by 
considering  these  answers  to  the  following  questions. 

(i)  What  are  we  studying?  The  free  energy  f  the  Lenz-Ising 
model  will  be  examined.  <2)  Why  are  we  stu  .ying  this  model? 
We  are  interested  in  obtaining  a  method  that  will  permit  the 
calculation  of  the  free  energy  for  systems  of  interacting 
particles.  The  Lenz-Ising  model  is  the  simplest  of  all  the 
models  of  interacting  particles.  No  exact  solutions  in 
three  dimensions  for  arbitrary  T,B  yet  exist.  (3)  How  shall 
we  study  this  model?  The  Morita  Cluster  Variation  Method 
(CVM)  will  be  used  to  calculate  an  approximate  value  for  the 
free  energy  of  this  model .  <  4 )  Why  do  we  use  this  method? 

Variational  methods  have  a  track  record  as  a  superior 
approximation  method.  The  Morita  Variational  Method  has 
been  found  to  behave  strangely;  we  wish  to  resolve  this 
strangeness . 

In  addition  to  obtaining  answers  to  the  above 
questions,  this  work  will  assist  in  gaining  insight  into  the 
phase  transitions  exhibited  by  various  systems  at  a  critical 
temperature  peculiar  to  each  system.  Examples  of  phase 
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changes  are  the  transition  from  a  non-ferromagnetic  to  a 
ferromagnetic  state  in  Fe  and  Ni  i.e.  the  long  range 
ordering  that  causes  the  material  to  become  magnetic.  Other 
examples  are  a  lattice  gas  and  a  binary  alloy  such  as  Cu-Zn. 

Application  of  the  CVM  to  Calculate 

the  Free  Energy  of  the  Lenz-Ising  Model 

One  way  to  study  the  behavior  of  an  N-particle 
system  is  to  analyze  this  behavior  into  the  behavior  of  the 
N  1-particle  subsystems,  the  N*(N-l)/2  2-particle 
subsystems,  the  N*(N-l)*(N-2)/(3! )  3-particle  subsystems, 
and  so  on.  This  method  replaces  the  study  of  one  object 
— the  N-particle  system — by  the  study  of  2N~1  objects.  Why 
would  one  do  this? 

One  always  does  this  for  systems  whose  particles 
(1)  do  not  interact  with  each  other  whence,  the  system’s 
behaviors  are  entirely  determined  by  the  N  1 -particle 
properties  and  (2)  are  all  similar  in  their  properties, 
whence,  all  N  1-particle  properties  are  the  same.  Then  this 
approach  replaces  the  study  of  an  N-particle  system  by  a 
1-particle  system.  Examples  of  this  approach  include  the 
standard  studies  of  the  ideal  classical  and  quantal  gases, 
and  the  ideal  paramagnet. 

This  approach  has  been  extended  to  include  those 
systems  of  interacting  particles  for  which  "normal 
coordinates"  have  been  located,  coordinates  in  which  there 
are  no  interactions.  The  main  example  here  is  the  set  of 
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vibrations  of  lattice  ions  in  the  harmonic  approximation, 
which  are  replaced  by  non-interacting  "quasiparticles”.  In 
this  case,  however,  the  1-quasiparticle  properties  are  not 
all  the  same;  each  class  must  be  treated  individually,  and 
the  final  results  obtained  by  summation. 

However,  arguments  that  "only  small  subsystems  are 
important"  are  rarely  rigorous  for  non-ideal  systems.  Nor 
do  all  of  these  sorts  of  arguments  include  estimates  of  the 
errors  made  in  ignoring  the  remaining  (enormously  many) 
subsystems.  Further,  it  is  often  difficult  to  see  how  to 
improve  a  treatment  of  this  type,  which  might  be  promising, 
but  not  quite  good  enough.  These  difficulties  all  come  from 
the  same  sources: 

1.  the  choice  of  which  subsystems  are  included  in  the 
study  is  usually  made  on  an  intuitive  basis,  and  not  on  the 
basis  of  a  general  principle,  and 

2.  ignored  subsystems  are  not  explicitly  considered  at 
all,  so  that  it  is  difficult  to  judge  the  error  that 
ignoring  them  incurs,  or  how  to  include  them  in  a  better 
approximation . 

The  usual  criterion  for  including  subsystems  is  that 
their  particles  are  all  physically  close,  and  in  intense 
interactions;  subsytems  whose  particles  are  all  distant,  and 
not  interacting,  are  usually  ignored.  This  is  why  the 
systems  studied  are  often  called  "clusters";  i.e.  a  number 
of  similar  things  in  close  proximity,  somewhat  separate  from 
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other  things.  There  does  not  seem  to  be  a  common  name  for 
particles  that  are  not  parts  of  clusters.  That  is,  common 
language  reflects  the  procedure  sketched  above:  we  know  how 
to  treat  N-particle  systems  that  obviously  separate  into 
essentially  independent  subsystems  (clusters),  but  we  have 
difficulty  when  this  separation  is  not  perfect. 

T.  Morita[l]  has  presented  a  formalism  for  the 
equilibrium  statistical  mechanical  study  of  N-particle 
systems  in  which  all  2N  subsytems  explicitly  appear. 

One  then  makes  an  explicit  choice  as  to  which  will  be 
retained  and  which  will  be  ignored  in  further  study.  This 
leads  to  an  approximate  expression  for  the  free  energy  of 
the  N-particle  system  in  terms  of  quantities  related  to  the 
retained  subsytems [ 2-4 ] .  Minimizing  this  approximate  free 
energy  by  varying  these  quantities: 

1.  gives  an  estimate  of  the  exact  equilibrium  free 
energy  ,  and 

2.  gives  estimates  for  all  equilibrium  thermal 
behaviors  of  the  system. 

It  had  been  thought  that  it  gave  an  upper  bound  to 
the  free  energy  --  we  now  know  that  this  is  wrong  1 

There  is  only  one  sort  of  approximation  involved  in 
this  procedure:  the  choice  of  which  subsystems  to  retain  for 
further  consideration,  and  which  to  ignore.  Many  other 
strategies  require  a  series  of  nested  approximations;  it  is 
frequently  difficult  to  determine  their  relative  affects  on 
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the  accuracy  of  the  result.  We  now  know  that  Morita’s 
method  is  also  of  this  uncontrolled  type! 

Morita’s  method  can  be  systematically  improved  by 
retaining  previously  ignored  subsystems.  Since  these 
explicitly  enter  the  original  expression  for  the  free 
energy,  it  is  possible  to  form  some  estimate  of  their 
importance.  It  is  easy  to  incorporate  any  "physical 
intuitions"  one  may  have  as  to  the  importance  of  certain 
clusters . 

Morita  has  called  his  procedure  the  CLUSTER  VARIATION 
METHOD.  We  shall  present  it  by  defining  a  "cluster"  in  a 
general  way,  one  that  does  not  rely  on  some  particular 
attribute  such  as  physical  proximity.  We  shall  restrict 
ourselves  to  classical  systems,  and  we  shall  suppose  that 
the  particles  are  arranged  on  a  crystal  lattice.  Then  we 
shall  define  the  microstates  and  the  macrostates  of 
clusters,  and  then  distinquish  between  cluster  functions 
and  cluster  quantities.  We  shall  introduce  the 
distinction  between  extrinsic  and  intrinsic  variables, 
recall  the  usual  development  of  the  Hamiltonian  into  a 
series  of  intrinsic  terms,  and  exhibit  Morita’s  novel 
development  of  the  entropy  into  such  a  series.  We  exhibit 
the  free  energy  of  the  system  as  a  sum  of  a  series  of 
intrinsic  terms,  each  due  to  a  particular  cluster.  This 
series  is  exact,  but,  as  it  contains  2N  -1  terms,  it  is 
too  difficult  for  direct  calculation  of  results  for 
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macroscopic  systems.  But  this  series  is  ideal  as  a  starting 
point  for  approximations  --  one  partitions  the  set  of 
clusters  into  those  which  one  will  retain  in  subsequent 
work,  and  those  which  one  will  ignore,  and  then  deletes  from 
the  sum  all  the  terms  relating  to  ignored  clusters.  In  this 
work  the  truncated  expansions  of  the  free  energy  are  also 
called  approximations.  Thus  the  expansion  retaining  only 
the  first  term  is  called  the  first  approximation  and  the 
expansion  retaining  only  the  first  and  second  terms  is 
called  the  second  approximation  and  so  on.  In  this  work  the 
free  energy  will  be  calculated  for  the  first  through  the 
fourth  approximation  using  the  Morita  expansion. 

The  calculation  of  the  free  energy  using  these 
concepts  can  be  outlined  as  follows.  The  Helmholtz  free 
energy  of  the  Lenz-Ising  model [5]  is  expressed  as  a  function 
of  its  macrostate,  F  =  F(Po ) ,  for  the  first  through  the 
fourth  approximation  using  the  Morita  expansion.  For  each 
approximation,  that  macrostate  Po  for  which  the  free  energy 
is  stationary  is  found,  using  either  the  calculus  methods  or 
the  Simplex  algorithm[6] .  The  stationary  free  energy  is 
obtained  with  the  calculus  methods  by  setting  the  first 
derivative  of  F  to  zero'.  (<5'F/6‘Po  )po  =  0  and  determing  the 
mini-max  nature  of  the  stationary  state  Po  with  the  second 
derivative.  For  the  higher  approximations  the  Simplex 
algorithm  is  used  to  find  the  stationary  free  energy. 

The  method  can  be  summarized  as:  Write  the  free 
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energy  (using  "the  Hamiltonian)  for  the  Lenz-Ising  model 
using  the  Morita  cluster  expansion  and  minimize  it. 

Remembering  that  the  equilibrium  free  energy  is  the  one 
which  is  the  least  of  all  possible  free  energies  (F) ,  the 
minimization  will  give  two  results:  (1)  the  equilibrium 
macrostate  Po  ,  and  (2)  the  equilibrium  free  energyy  Fo  = 

F(Po).  This  method  is  completely  analogous  to  the  quantum 
mechanics  procedure:  minimizing  the  expectation  of  the 
energy  <E>  =  (psi,Hpsi)  gives  not  only  the  correct  ground 
state  energy  of  the  system  Eo  but  also  the  correct  ground 
state  psio . 

The  results  for  the  Lenz-Ising  model  in  the  presence 
of  a  magnetic  field  will  be  calculated  with  the  Morita  method 
and  compared  with  exact  calculations  for  systems  which  contain 
only  N=l,2,3,4  particles  (see  Chapter  III). 

Boltzmann.  Gibbs.  Morita 

Boltzmann  expressed  the  entropy  of  composite  systems 
whose  sub-systems  do  not  interact  as[7], 

Sb  =  -kBZ  i  =  i  N  tri  (  i )  Pi  (  i )  InPi  (  l ) 
where  kB  is  Boltzmann’s  constant 

tri  (  1 )  is  a  (generalized)  summation  over  all  the 
microstates  of  the  ith  sub-system, 

Pi  < 1 )  is -the  probability  law  for  the  microstates  of 
the  ith  sub-system,  that  is,  Pi  (  i )  is  the 
MACROSTATE  of  the  ith  SUB-SYSTEM. 

Later  this  was  extended  by  Gibbs[8]  to  calculate  the 
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entropy  for  composite  systems  with  arbitrarily  interacting 
subsystems  as , 

Sg  =  kBtr<  N)  P(  N>  inP<  N) 

where , 

tr< N)  is  a  (generalized)  sum  over  the  microstates 
of  the  composite  system, 

P<  N)  is  the  probability  law  for  the  macrostates 

of  the  composite  system,  that  is,  P(  N)  is  the 
MACROSTATE  of  the  COMPOSITE  SYSTEM. 

For  a  composite  system  of  many  interacting  particles 

the  calculation  of  Sg  is  prohibitively  difficult.  Morita 

expresses  the  total  or  extrinsic  entropy  of  the  composite 

system  as  the  sum  of  the  entropies  associated  with  each 

cluster  of  the  composite  system.  In  the  cluster  expansion 

he  refers  to  the  extrinsic  entropy  of  the  complete  system  as 

given  by, 

Sg  =  S(  N) extrinsic  =  S~< 1 > intrinsic  +  S~ (  2 )  intrinsic 

+  .  .  .S~<  N)  int. 

where,  is  the  intrinsic  entropy  of  all  particles 

taken  one  at  a  time  (singlets)  with  no  interaction  between 

the  particles  and  would  be  the  same  as  Boltzmann’s  Sb . 

Therefore,  Sb  =  S~(x) .  The  next  term  S~(2)  is  the 

intrinsic  entropy  of  a  pair  of  interacting  particles  and 

then  S~(3)  is  the  intrinsic  entropy  of  a  triple  of 

interacting  particles  and  so  on.  This  resolution  of  the 

entropy  also  clarifies  the  relation  of  Sg  and  Sb  since  we 

could  write  the  expansion  as 

S(  N )  =  Sg  =  Sb  +  S~  (  2  )  +  S~  (  2  )  +  S~  (  2  )  +  ...  S~  (  n  )  . 

« 

If  there  are  no  correlations  among  the  particles,  i.e. ,  if 
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P<  »)  =  Si  =lNP(  1  )  i  , 

then  S(n)  =0  for  each  n  >  1, 

and  Sg  =  Sb  . 

Then  Gibbs'  entropy  reduces  to  Boltzmann's  entropy.  If 
there  are  interactions  among  subsystems  (particles),  then  Sg 
will  differ  from  Sb . 

The  complete  Morita  expansion  for  the  free  energy 
(F)  is  exact  but  typically  we  only  use  the  first  few  terms. 

We  refer  to  the  terms  as  clusters  or  approximations  i.e.,  1st 
approximation  or  1-cluster  etc. 

Morita  decomposes  a  system  of  N  particles  into  the 
collection  of  all  subsets  of  these  particles: 

"1-clusters"  are  sets,  each  containing  1  particle;  there 
are  N  of  these. 

"2-clusters"  are  sets,  each  containing  2  particles;  there 
are  N(N-l)/2  of  these. 

"n-clusters"  are  sets,  each  containing  n  particles;  there 
are  N!/(N-n)!n!  of  these. 


The  entire  system,  then  is  "the  N-cluster" ;  there  is 
one  of  these.  We  could  use  the  familiar  notation  for 
"combinatorial  factor" : 

Cn(N)  =  N!/(N-n) !n! 

There  are  2N~1  clusters  in  all,  for  a  system  of  N  particles. 
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Any  particular  n-cluster  is  labeled  by  the  particles  it 
contains : 

C(n)il,i2,i3,.. .in.  The  microstate  of  the  n-cluster 

C(n)ii . in  is  Pni  1  ,  .  .  .  ,  i  n  =  (Mil, ...Pin).  The 

microstate  of  the  ith  particle  is  Pi  ;  the  microstate  of 
the  whole  system  is  PN  =  (Pi,P2,  ...  Pn).  The 

macrostate  of  the  whole  system  is  P<N)  =  probability  of 
each  value  of  pN.  The  macrostate  of  an  n-cluster  is 
P(  n )  i  l  ,  .  .  .  i  n  =  EP(  N)  where  sigma  is 
summed  over  all  the  microstates  of  the  particles  not 
in  the  n-cluster. 

The  following  approximations  from  Morita’s  expansion 
for  the  free  energy  are  the  expressions  that  we  will 
minimize.  Only  the  first  two  approximations  are  shown  here. 

First  Approximation  (first  term)  (see  also  p.41,68) 

This  first  term  in  the  Morita  cluster  expansion, 
known  as  the  first  approximation,  is  the  free  energy  for 
single  clusters-one  particle  at  a  time,  with  no  interaction 
between  clusters.  All  correlations  at  the  pair  level  and  up 
are  ignored. 

F  =  U  -  TS 
U  =  <H> 

H  =  -Po  (Si  =  iNXi  )B  -  y2J(Zi  =iNxi  2  )Z 
U  =  N(-EPoBx-J£ZJx2  ) 

where  F  =  free  energy  T  =  temperature 

U  =  internal  energy  S  =  entropy 
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Z  =  coordination  number*  B  =  external  magnetic  field 


H  =  Hamiltonian 
x  =  magnetization 
per  particle 

*Z  is  the  number  of  n.n. 
n.n.  =  nearst  neighbors 


J  =  exchange  integral 
JJO  =  magnetic  moment 
x2  =  interaction  between  n.n. 
(see  p.48) 

Z  =  2  in  one  dimension 
Z  =  4  in  2-dim  (square  lattice) 

Z  =  12  in  3-dim 

In  this  study  of  the  free  energy  using  the  Morita 
cluster  expansion  the  results  will  be  obtained  for  an  FGG 
lattice[9].  A  unit  cell  of  an  FCC  lattice  is  shown  in  Fig.l. 
The  Morita  expansion  for  only  the  first  term  is, 

S(  N)  =  S~( 1 ) 


where 

S~d>  =  SO)  =  2[all  singlets ]S<  i  >  i  . 
and  S(  1 )  i  =  kB  tr(  1 )  i  P(  1 )  i  lnF<  1 )  i  is  the  entropy  of 
each  single  subsystem  (or  particles)  and  S(  1  >  is  the  entropy 
of  all  the  N  subsystems,  each  taken  one  at  a  time,  so  that 
there  are  no  interactions  among  the  systems,  and  = 

S(D  is  the  INTRINSIC  entropy  of  these  NONINTERACTING 
subsystems . 

S  =  -kBN[Jg(l+x)ln(  .  )  +  Jg(l-x)ln(  .  )] 

F  =  N{-Po Bx-^ZJx2  +  kBT[J*(l+x)ln(  .  )+JS(l-x)ln(  .  )]} 
or  F/NkBT  =  (5F/N  =  -(PoB/kBT)x  -  (ZJx2/2kBT) 

+  K(l+x)ln(  .  )+>S(l-x)ln(  .  )] 

When  the  factor  that  precedes  the  natural  log  is  repeated  as 
its  function  then  the  notation  ln( . )  is  used. 
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Fig.  1.  THE  UNIT  CELL  FOR  A  FACE  CENTERED  CUBIC  (FCC) 
LATTICE.  THE  NUMBER  OF  NEAREST  NEIGHBORS  IS  12. 
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§  =  GF/N  =  -(pofJB)x  -  (5^,JZ)x2  +  [J$(l+x)ln(  .  )+J$(l-x)ln(  .  )] 
Using  a  reduced  "temperature  Tr  an  alternate  form  can  be 
obtained  ((i=l/kBT) 

Tr  =  kBT/J  =  l/fl  J 

Then  P o B/kB T  =  PoB/kBT*J/J  =  PoB/JTr  =  Br /Tr 
where  Br  =  P  o  B/ J 

5  =  -(Br/Tr)x  -  (Z/2Tr)x2  +  J$(l+x)ln(.)  +  J$(l-x)ln(.) 

For  Br  =0 

5  =  -(Z/2Tr)x2  +  %( 1+x)  ln(  .  )  +  *(l-x)ln(.) 

We  find  the  "best"  values  for  x  by  minimizing  5 . 

0  =  d5/dx  =  - ( Z/Tr ) x  +  %ln%(l+x)  -  %ln%(l-x) 

=  - (  Z/Tr  ) x  +  3£ln[(l+x)/(l-x)] 
x  =  %(Tr/Z)ln[(l+x)/(l-x)] 

Second  Approximation  (second  term)  or  2-clusters 

+  1-cluster 

The  second  approximation  is  the  sum  of  the  first  term 
(first  approximation)  plus  the  2-cluster  term  of 
Morita’ s  expansion, 

S(  « )  =  S~  <  i )  +  S~  (  2  ) 

where 

S~(2>  =  I  (all  pairs  )  S~  <  2  )  i  ,  j 

and  S~(2)i,j  =  S(2)i,j  -  S~(Di  -  S~d)j  is  the  INTRINSIC 
entropy  for  subsystems  of  pairs  of  particles  S~d)i  and 
S~(i)j,  including  any  interactions  they  may  have,  and 
S(  2  >  i  ,  j  =  -kB  tr<  2  )  i  ,  j  P(  2 )  i  ,  j  lnP(  2 )  i  t  j  is  the  entropy 
of  S~(!)i  and  S~d)j  including  interactions.  The  entropy 
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due  only  to  the  single  particles  S"(1)  is  subtracted  from 
the  pair  entropy. 

F  =  U  -  TS 

U  =  -(NMoB)x  -  (^NZJ)y 

S  =  S~  (  1 )  +  S'-  (  2  ) 

=  -kBN{^(l+x)ln^(l+x)  +  ^(l-x)ln^(l-x) 

+  ^Z[^(l  +  2x+y)ln34(l  +  2x+y)  +  2**(  l-y)ln*(l-y) 

+%( l-2x+y ) ln%(  l-2x+y )  ] 

-2[^(l+x)ln^(l+x)  +  ^(l-x)ln^(l-x)]}. 

F  =  N{-PoBx  -  ^ZJy  +  TkB  [^(l±x)  ln^(l±x) 

+  ^Z{^(l±2x+y)ln54(l±2x+y)  +  %(  1-y )  ln*4(  1-y ) 

-  d±x)lnJ£d±x)}]}. 

F/NkBT  =  5  =  -(Po^B)x  -  ^JJZy  +  ^(l±x)ln^(l±x) 

+  J£Z{%(  l±2x+y )  ln(  .  )  +  ^(l-y)lnH(l-y) 

-  ( l±x)  ln%(  l±x)  }  . 

5  =  -  ( Br  /Tr  )  x  -  ^(Z/Tr)y  +  ^(l±x)ln^(l±x) 

+  ^Z{^(l±2x+y)ln(. )  +  ^(l-y)ln^(l-y) 

-  (l±x)ln^(l±x)}. 

For  Br  =  0 

3  =  -^(Z/Tr)y  +  %(l±x)lnH(l±x)  +  ^Z[?4d±2x+y ) ln(  .  ) 

+  %( l-y )  ln^(  1-y )  -  ( l±x)  ln^(  l±x)  ] 

0  =  fi  (x,y)  =  d3/dx  ;  0  =  f2(x,y)  =  dS  /dy 

It  is  helpful  to  recall  that  the  first  few  terms  of 
the  Morita  expansion  are  known  by  other  names  -  the  names  of 
their  inventors,  and  later  reappeared  as  terms  in  the 
cluster  expansion  for  the  free  energy.  (See  the  names  of  the 
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inventors  associated  with  the  approximations  shown  in 
figure  5,  Chapter  II.)  This  is  interesting  since  in  their 
original  form  it  was  sometimes  difficult  to  see  how  their 
inventors  were  able  to  logically  derive  the  equations  for 
the  free  energy.  As  these  results  are  rederived  in  the 
Morita  expansion  they  appear  as  a  consequence  of  a  very 
orderly  and  logical  development  and  the  relation  of  the 
approximations  to  each  other  is  now  understood.  The  physics 
is  clear  in  this  new  derivation  as  compared  to  the  original 
derivations  where  it  was  sometimes  obscure.  The  first  term 
or  first  approximation,  known  as  the  1-cluster  term  in  the 
Morita  expansion,  is  the  same  as  the  Weiss  model  which  uses 
the  self-consistent  field  or  molecular  field  approximation 
for  f erromagnets .  The  Bragg-Williams  model  (1934)  is  a  more 
rigorous  statement  of  the  Weiss  model  and  gives  similar 
results,  (see  figure  5).  The  second  term  is  the  2-cluster 
term  and  was  originally  known  as  the  Bethe-Pierls 
approximation  used  to  describe  Cu-Zn  alloys  and  is  an 
improvement  over  the  Bragg-Williams  approximation.  The 
Bethe-Pierls  model  (1936)  takes  into  account  specific  short- 
range  order,  i.e.  local  correlations  between  spins.  The 
second  approximation  is  the  sum  of  the  1-cluster  and  2- 
cluster  terms. 

Some  items  of  procedural  interest 

A.  Particle  Names 

Label  the  N  particles  of  a  system  using  the  N 
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counting:  1,2,...,N.  This  labeling  is  arbitrary.  The  same 


same  labeling  is  used  unchanged  throughout  this  work. 

A  "cluster"  of  particles  is  a  set  of  particles, 
that  is,  an  unordered  list  of  particles.  For  example, 
{il, .  .  .  ,  in}  =  {i«}na=i  is  a  particular  cluster  of  n 
particles.  This  is  also  called  an  n-cluster  of  the 
particles  of  the  system.  The  individual  particles  in  that 
cluster  are  called  il  ,  i2  ,  .  .  .  in  .  In  a  system  of  N 
particles,  there  are  the  following  n-clusters: 
size  of  cluster  generic  name  specific  names 

N  1-clusters  {i}  {1} , {2} , . . . , {N} 

^N(N-l)  2-clusters  {il,i2>  {1,2} 

(1/6)N(N-1) (N-2)  3-clusters  {ii,i2,i3}  ... 


N(N-1 )  N-cluster  {il , . . . , in }  ... 

In  all,  there  are  2N  -1  n-clusters. 

Cluster  function: 

A  "cluster  function"  is  a  function  of  the  microstate 
of  the  particles  making  up  a  cluster.  For  example, 
f(N)(il,...,in)  =  f  (Pi  l  ,  Pi  2  ,  Pi  3  ,  .  .  .  ,  Pi  n  )  is  an 

n-cluster  function,  where  Pi  is  the  microstate  of  the  ith 
particle  and  Pi  1  ,Pi 2  ,  .  .  .  ,Pi n  is  the  microstate  of 
{il  ,  .  .  .  ,  in  }  . 

So  Pi  is  the  microstate  of  particle  ii 
P  2  "  "  "  "  "  i2 
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In  particular,  the  macrostate  of  an  n-cluster  is  a  cluster 
function. 

Names : 

There  are  N  particles  arranged  on  N  sites  of  a 
lattice.  Sites  are  located  at  the  positions 
r  =  ni  bi  +  n2  b2  +  ...  n3  b3 

where 

bi  is  the  itl>  base  vector  of  a  unit  cell  of 
the  lattice,  and  each  ni  is  an  integer. 

Sometimes  it  will  be  convenient  to  name  a  particle 
for  the  site  it  occupies.  Then  we  speak  of  the  particle 
ni,n2.  However,  for  much  of  the  formal  discussion  of  this 
work,  it  will  be  convenient  to  assign  the  cardinal  numbers 
{1,2,...,N}  to  the  particles  in  some  way,  and  then  denote  a 
particle  by  this  cardinal  number.  So,  we  speak  of  the  it h 
particle . 

P<n>il,i2,...,in  =  P(  n)  i  1  ,  i  2,  .  .  .  in  (Pi  1  ,  Mi  2,  ...Pin) 
is  the  macrostate  of  the  cluster  C(n)ii,  ...in,i.e., 

P(  »)ii,  . . . ,in  is  the  probability  of  the  microstate 
Mil . Pin. 

A  "cluster  quantity"  is  the  average  of  a  cluster 
function;  it  depends  on  the  macrostate  of  the  cluster. 

A  "cluster  variable"  is  either  a  "cluster 
function" ,  or  a  "cluster  quantity" .  Cluster  variables 
can  be  "extrinsic"  or  "intrinsic".  We  suppose  we  can  obtain, 
perhaps  by  a  direct  measurement,  the  value  of  a  cluster 
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variable  for  any  isolated  cluster:  these  are  extrinsic, 
then  define  the  intrinsic  cluster  variables: 
f(  i)  i  =  f*< 1)  i 

f<2)id  =  (f*(Di  +  r<i)j  )  +  f‘(2>ij 

f(3)ijk  =  (f‘(i)i  +  f*(i)j  +  f"(Dk  ) 

+  (f~<  2>  i  j  +  f*(  2)  ik  +  f*(  2)  jk  ) 
+(f‘(3)ijk  ) 


f<  n>  i  X  ,  i  2,  .  .  .  in  =  (E  j=inf*(  1  )  i  j  ) 

+  (Sjkf-<  2)  i  jik  ) 


These  can  be  inverted  to  define  any  particular  intrinsic 
cluster  variable  entirely  in  terms  of  extrinsic  cluster 
variables  only. 

f*(  Di  =  fd)i 

f“<  2)i  i  =  f<  2)  i  j  -(f<  Di  +f<  l )  j  ) 

f‘(3)ijk  =  f  (  3  )  i  j  k  -  (f<  2)  i  j  +f(  2)  ik+f(  2)  jk  ) 

+  (f(  Di  +f(  1  >  j+f<  1  >k  ) 


We 
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f  (n)il,i2,...in  -  f  (  n  )  i  1  ,  i  2  ,  .  .  .  i  a  ~  (J-  f(  n  1  )  i  .  .  .  i  ) 

+  (Ef<  n-2)  i.  .  .  i  )  -  ... 

Examples  : 

H  =  EhA(  1  )  i  +  Eh~  <  2>  i  j  +...Eh<N)i2...N 

The  Hamiltonian  H  is  an  extrinsic,  N-cluster  function. 

It  is  expanded  into  intrinsic  n-cluster  functions. 

U  =  E  <hA(  1  )  i  >  +  £  <h*<  2)  i  j  >  +  .  .  . 

The  internal  energy  U  is  an  extrinsic  N-cluster  quantity. 
It  is  expanded  into  intrinsic  n-cluster  quantities. 

S  =  -kBZMlZP2  .  .  .2PnP<  N)  lnP(  H> 

=  S*<  1 )  +  S~<  2)  +  ...  +  S~<  **) 
where 

S~(  n  =  EiSi~<  i) 

SA(  2)  :  HijSA(2)ij 

S  (3)  z  £  i  j  k  S  (  ^  )  i  j  k 

The  entropy  S  is  an  extrinsic  N-cluster  quantity.  It  is 
expanded  into 

S(Di  =  SMD  z  -kBl  (Pi  )P(  i  )  i  lnP(  i  )  i 
S(  2)  i  j  =  SC  1 )  i  +  S(  1 )  4  +  S~<  2)  i  4 
=  -kBlPiZP2P(  2)i  jlnPC  2)1  j 

etc . 

B.  Summing  over  clusters,  and  the  effects  of  symmetries 

An  n-cluster  is  defined  by  the  n  particles  it 
contains:  the  order  of  these  particles  is  of  no  consequence 
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to  "the  definition  of  the  n-cluster,  [So  {1,2}  is  the  same  as 
{2,1},  and  { lot  } (oc  =  l-->n) ,  is  the  same  cluster  as 
{Pia  ) }  (oc  =  l  —  >n) ,  where  P  is  any  one  of  the  ni 
permutations  of  the  n  numbers  1,2, ... ,n  ]  So,  while  there 
are  n!  distinct  lists,  i.e.  arrangements,  of  the  n  numbers 
l,2,...,n,  all  these  lists  name  the  same  cluster'- 
{1,2, ...,n},  { io<  }oc . 

Sums  of  1-cluster  quantities: 

E  V(  i  }  E  (  a  1  1  1-clu.sters)  Q(l)i 

=  Zi  =  lNQ(  1  )  i 

=Q(  1 )  l  +  Q<  i  >  2  + .  .  .  +  Q<  i )  3 

Z  Vi  j  Zall  2-cluateraQ(2)ij 

=  ^Zi=l«Z  j=lNQ<  2)  i  j 

(i  =/  0) 

=  J$(Q<  2)12  +.  .  .+  Q<  2)  in 
Q(  2  )  2  1  +  .  .  .  +  QN 1 


Q(  2  )  n  1  +Q<2)N2  +...  ) 

=  Zi=lN-lE  j=i+lNQ(  2)  i  j 
=  Q(  2  )  i  2  +  Q(  2  )  i  3  +...  +  Q(  2  )  i  n 

+  .  .+  QN<  2)  N_i 

[Note  that  Q(2)ij  =  Q(2)ji] 

Zvi  jk 

Zall  3-clusters  Q(2)ijic 
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=•(1/3!  )IiS  jXkQ<  3)  i  jk 
(i=/j=/k) 

r£i=lN-2Ej=i+iN-lEk=j+lNQ(3)ijk 

[Note  that  Q<3>ijk  =  Q(3)ikj  =  Q(  3 )  j i  k  =  Q(3>jki 

=  Q(  3)  ki  j  =  Q(  3)  kji  ] 


Sail  n-clustersQ(n)il,i2...in 

=  (  1/n  !  )Eil=lNEi2  =  lN  .  .  . 

E  i  n  =  1  N  Q(  n)  i  1  i  2.  .  .  in 

( ia  /=  id  for  any  <x  ,d  -  1,N  except  when  <x  =  d 

— L  il=l^"nEi2  =  il+lN“n  +  l 

Iin  =  in-1+1^Q(  n  )  i  1  i  2  .  .  .  in 


For  similar  particles  on  a  regular  lattice,  we  can 
drastically  simplify  the  above  sums.  We  use  the 
"n-coordination  numbers'1,  defined  as  follows: 

Z( 2 ) p  =  pair  coordination  number  for  spacing  p 

=  number  of  2-clusters  whose  spacing  is  p, 

Z(  3 ) pqr  =  tri-coordination  number  for  spacings  p,q,r 
=  number  of  3-clusters  whose  sides  have 
lengths  p,q,r,  with  piqir. 

Z(  4) pqrs  =  quad-coordination  number.  .  . 
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etc . 


Tall  1  -  el  uster  sQ(  ^  )  i  —  N*Q(  1  )  <X 

for  a  =  any  one  of  {1,2, ... ,N}  . 

Sail  2-elu3tersQ(2)ij  —  p  Z(  2  )  p  Q<  2  )  (x  ,  (1  , 

for  a  ,(i  ,  -any  of  {1 , 2,  .  .  .  N}  such  that 
1  r_>*  -  r_>(l  !  — >p 

Sail  3-clusteraQ(3)ijit  —  (1/3!)  NE  pZqZr  Z(  2)  p  q  r  Q(  ®  )  QC ,  (i  ,  P 

for  any  «,(3,r  or  any  permutation  of  these, 
etc . 

C . Microstates  and  Macrostates 

The  microstate  of  the  jth  particle  is  Pi  :  its  values 
fall  in  some  discrete  or  continuous  set.  For  the  Lenz-Ising 
model,  the  values  are  +1  and  -1.  The  microstate  of  an 
n-cluster  is  the  set  of  the  microstates  of  each  particle 
contained  in  it-- 

M<n>il,...in  =  {Pi  1  ,  Pi  2  ,  .  .  .Pi  n  } 

In  particular,  the  microstate  of  the  N-particle  system  is 
the  set  of  the  microstates  of  each  particle: 

p<  n)  =  {pi  , p 2  ,  .  .  .  ,Pn  }  . 

The  macrostate  of  the  N-particle  system  is  a 
probability  law  defined  on  the  set  of  possible  microstates 

p<  H)  : 

P(N>  (P<N)  )  =  probability  of  P<  N> 

The  task  of  equilibrium  statistical  mechanics  is  to 
determine  this  macrostate. 

From  the  macrostate  of  the  system,  P(  N)  ,  we  can  extract 
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the  macrostate  of  any  subsystem: 

P<  n>  (Mi  1  ,  .  .  .  ,Pi  n  ) 

=  SPi...EPNP(N)(Pl,P2,...,PN) 

( except  Pi  ,  P  2  ,  .  .  .Pin) 

(1)  From  this,  it  follows  that  the  macrostate  of  any 
n-cluster  determines  the  macrostate  of  all  the  n' -clusters 
(  n’ <  n  )  contained  within  it. 

(2)  We  note  that  all  macrostates  are  normalized: 

EPil...EPinP(n)il . in  =  1 

(for  all  n  and  all  {ia}[a  =  l  ->  n]  )  . 

D.  A  note  on  summing  over  n-clusters 

We  often  sum  over  all  the  various  n-clusters,  for  a 
given  value  of  n.  We  also  often  sum  over  all  the  various 
microstates  of  a  given  n-cluster  (e.g.,  when  computing  a 
thermal  average).  It  is  convenient  to  represent  these  very 
different  sorts  of  sums  with  very  different  sorts  of 
notations.  So,  henceforth  we  will  use 

Si lEi 2  .  .  .Si n  [  . ]  =  Sail  n-cluater  [.] 

( ia  =  /  ifi  Vtx ,  (3  xxx  a  =(5  ) 

SPilSPi2  .  .  . SPi  n  [  .  .]  =  tr(  “)  i  1 ,  i  2.  .  .  in  [  .  ] 
respectively,  for  the  two  different  sorts  of  sums. 

E .  Representation  of  macrostates 

It  is  a  well-known  result  that  a  probability  law  is 
equivalent  with  the  set  of  all  of  its  moments. 

A  common  (sketch  of  a)  proof  is: 
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(1) 


P(  x )  =  probability  law  for  a  random 


variable  x. 


(2) 


R(k)  =  fl-®®P(x)exp(ikx)dx 


is  the  Fourier  transform  of  P(x) . 


[H  represents  an  integral  sign] 


(3) 


R(k)  =  R(  0 )  +  R’(0)k  +  R”(0)k2/2  +  ... 


is  the  Taylor’s  series  of  R(k)  . 
(4)  By  inspection  of  (2)  : 

R( 0  )  =  f1-®®P(x)dx  =1 

by  normalization  of  P(x). 

R’(0)  =  if1-®®xP(x)dx  =  x-  is 
the  l3  moment  of  x . 

Note:  x“  represents  x  with  a  "bar"  superscript. 
R’’(0)  =  (  i2  )fl-®®x2  P(x)dx 

=  (X2)- 

is  the  2nd  moment  of  x. 


(5)  Hence,  R(k)  =  1  +  ikx-  +  ^(i2)k2x2"  +  ... 

(6)  Hence,  if  we  are  given  P(x),  then  we  can  calculate 


all  the  moments,  either  directly, 


(x^)-  =  fl-®®x^P(x)dx, 


or  indirectly, 


(xn)~  =  (dnR/dkn)k  =  o 


If  we  are  given  all  the  moments,  then  we  can  construct 
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R(k)  =  1  +  ikx-  +  (i2  )k2  (x2  )-  +  ...  , 

and  inverse-Fourier  transform  it  to  get  P(x) 


P(x)  =  ^rrn-®®exp( -ikx)R(k)dk 


QED. 


This  rule  extends  to  the  case  of  a  set  of  (possible) 


correlated  random  variables  : 


P(x)  is  equivalent  with  na=inxia  for  all  n. 


When  the  microstate  of  each  particle  is  bi-valued, 
Mi  =  ±1,  then  the  macrostate  of  the  particle  is 


where 


x<Di  =  <Mi  >  =  tr<  i  )  i  [P<  i )  i  (Mi  )Mi  ]  . 

clearly,  P(l)i(Mi)  is  equivalent  with  <M2>  each 
determines  the  other. 

The  macrostate  of  larger  clusters  can  also  be  expressed  in 
terms  of  moments : 


1  +  X<  1  )  i  -  x<  1 )  j  +  x<  2  >  i  j 


X  T  X\  Ml  -  XV  1  )  J  T  xi  6  /  1  J 

1  -  X<  1  )  i  +  X<  1  )  j  -  X<  2  )  i  j 

1  -  X<l)i  -  x(1)j  +  X(  2  )  i  j 


where 


x<2)ij  =  <M  i  M  j  >  =  tr(  2  >  i  j  [P<  2  )  i  j  (Mi  ,M  j  )Mi  M  j  ] 


P(  3)  i  jk  (MiM  jMk  )  =  P(  3)  i  jk  1  1  1 


11-1 
1-11 
1  -1  -1 

-111 
-1  1  -1 

-1  -1  1 

-1  -1  -1 
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=  1/8 


l+x<  i ) 
1+ 

1  + 

1+ 

1- 

1- 

1- 

1- 


i  +X< 
+ 


+ 

+ 


1  ) 


+X<  1  )  k  +  X<  2  )  i  j  +x<  2)  ik  +x<  2  )  j  k  +x<  3  )  i  j  k 
+  -  + 


+ 

+  + 
+ 


+ 

+ 


+ 

+ 


+ 


+ 

+ 

+ 


where 

X<  3)  i  jit  =  <MiJJjPk> 

=  tr<  3)  i  jk  [P<  3  >  i  jit  (pAp  jpk  )pip  jpk  ] 

And  so  on . 
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CHAPTER  II 


SHORT  HISTORY  OF  THE  LENZ-ISING  MODEL 

The  following  remarks  are  intended  to  provide  some 
additional  background  information  for  the  Lenz-Ising  model 
with  particular  reference  to  ferromagnetic  systems. 

We  focus  attention  on  the  Lenz-Ising  model  as  a 
mathematical  object  existing  independently  of  any  particular 
physical  system.  The  Hamiltonian  used  in  the  Lenz-Ising 
model  is  a  definition  of  this  model.  The  Hamiltonian  we  are 
using  is  for  nearest  neighbor  (n.n.)  interaction  only 
between  the  particles  in  the  system.  This  system  for 
interacting  particles,  is  the  simplest  one  that  it  is 
possible  to  study.  It  is  not  an  exact  description  of  any 
specific  physical  system  such  as  a  ferromagnet,  an  alloy  or 
a  liquid.  The  results  can  be  expected  to  give  the  general 
shape  of  the  magnetization  curve  and  the  value  of  the 
critical  temperature  but  not  the  details  for  any  particular 
physical  system.  This  however  is  only  a  limitation  of  the 
Hamiltonian  that  is  used  and  the  mathematics  necessary  to 
study  it.  Results  that  more  accurately  describe  a  specific 
physical  system  can  be  obtained  by  using  a  Hamiltonian  that 
is  an  improved  representation  of  that  system. 
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History 


To  explain  the  experimental  results  obtained  by 
Stern  and  Gerlach  in  1922  it  was  proposed  by  Uhlenbeck  and 
Goudsmit  in  1925  that  the  electron  possesses  a  spin  S=l/2 
and  that  it  is  oriented  either  parallel  or  anti-parallel  to 
the  f ield[ 10 ] .  Ising  published  in  the  same  year  the  results 

for  a  model  based  on  a  suggestion  of  his  thesis  advisor 
Wilhelm  Lenz,  that  if  electrons  were  located  on  a  lattice 
and  if  an  interaction  were  introduced  between  nearest 
neighbor  spins  that  favored  parallel  alignment  of  spins , 
then  at  sufficiently  low  temperatures  the  spins  would  all  be 
parallel,  and  the  model  might  provide  an  atomic  description 
of  ferromagnetism [11 , 12] .  The  Hamiltonian  corresponding  to 
this  model,  and  which  expresses  the  internal  energy  of  the 
system  can  be  written  in  the  form, 

H  =  -PoBEPi  -  JEij(sum  over  all  n.n.  pairs  )PiPj 
where  the  interaction  integral  (J  >  0)  represents  the 
interaction  between  spins,  B  is  the  magnetic  field,  Po  the 
magnetic  moment  of  a  single  spin,  and  Pi  or  Pj  =±lis 
the  spin  state  of  the  ith  electron  on  the  lattice.  The 
suffix  i  runs  over  all  sites  of  the  lattice,  and  <i,j> 
over  all  pairs  of  sites  i  and  j  which  are  nearest  neighbors. 

Ising  studied  this  model  in  one  dimension;  his  exact 
results  showed  that  spontaneous  magnetization  did  not  occur 
at  any  temperature.  This  result  is  correct.  He  extended 
the  one  dimensional  result's  to  higher  dimensions, 
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incorrectly,  and  did  not  find  any  phase  change  in  two  or 
three  dimensions.  This  result  is  incorrect.  Subsequent 
calculations  by  others  has  shown  that  there  is  not  a  phase 
transition  in  one  dimension  as  correctly  shown  by  Ising,  but 
for  two  or  more  dimensions,  spontaneous  magnetization  does 
occur  below  a  critical  temperature.  It  was  unfortunate  that 
Ising’ s  results  in  two  and  three  dimensions  were  incorrect 
since  it  caused  the  model  to  be  discarded  for  many  years. 

The  model  was  later  rediscovered  and  solved  exactly,  in  zero 
magnetic  field,  for  two  dimensions  by  Onsager  in  1944[13,14], 

Lenz  -  Ising  Model 

Further  work  has  shown  that  by  a  change  of  names  the 
Lenz-Ising  model  can  be  made  to  simulate  various  systems 
[15,16].  Examples  of  such  systems  are  (1)  magnets,  in 
which  each  molecule  has  a  "spin"  that  can  be  oriented  up  or 
down  relative  to  the  direction  of  an  external  applied  field; 
(2)  binary  alloys  such  as  Cu-Sn;  (3)  liquids  which  can  be 
represented  by  molecules  and  "holes"  (i.e.  empty  spaces)  on 
a  lattice  (this  is  called  a  "lattice  fluid")  and  there  are 
more  examples.  These  physical  systems  can  all  be 
represented  abstractly  by  the  same  model  which  can  be 
described  in  the  following  way.  We  assign  a  two-valued 
variable  (+1  or  -1)  to  each  node  of  a  regular  space  lattice 
to  represent  the  spin  Mi  ,  associated  with  each  node  i  of  the 
lattice.  For  a  magnet  M  corresponds  to  an  electron  spin 
state.  For  an  alloy  M  corresponds  to  an  ion  type  (say  Cu 
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vs.  Sn)  .  For  a  "lattice  fluid"  JJ  corresponds  to  the 
presence  or  absence  of  an  atom.  A  configuration  of  the 
lattice  is  a  particular  set  of  values  of  all  the  spins,  and 
for  N  nodes  there  are  2N  different  configurations.  An 
example  of  a  configuration  is  shown  in  Fig. 2. 

In  this  model  it  is  assumed  that  the  forces  between 
molecules  are  only  short  range  i.e.,  only  nearest  neighbor 
(n.n.).  When  neighboring  spins  are  the  same  (both  +1  or 
both  -1)  the  energy  is  -J  and  when  they  are  different  (one 
is  +1,  the  other  is  -1)  the  energy  is  +J[12] . 

The  interaction  tends  to  align  the  n.n.  spins  as 
parallel  and  to  give  aligned  spins  the  lower  energy  state. 
Heisenberg  in  1926  was  the  first  one  to  describe  the  reason 
for  this  interaction  between  spins.  It  is  an  electrostatic 
interaction,  a  Coulombic  type  force  that  causes  the  spins  to 
align.  And  so  the  spin  couplings  which  manifest  themselves 
as  magnetic  effects  are  in  fact  due  to  electrical  forces 
which  cause  the  spins  to  be  parallel  and  antiparallel.  The 
symmetric  or  anti -symmetric  wave  functions  that  the  electron 
spins  are  part  of  must  of  course  satisfy  the  rules  of 
quantum  mechanics.  For  the  three  physical  systems  described 
above,  this  interaction  could  result  in  (1)  spontaneous 
magnetization,  with  most  spins  in  the  same  direction,  even 
with  B  =  0;  (2)  transition  between  an  ordered  superlattice 

and  a  lattice  with  random  arrangement  of  atoms  on  the 
lattice  points;  (3)  condensation  of  molecules  in  one  region 
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THE  ENERGY  OF  THE  CONFIGURATION  IS, 

H  =  ~mb  B  I  /ij  -J I  jij  jij 
i  U 

(n.n.  PAIRS) 

H  =  -mb  B  (+6-3)  +  (-7J  +  5J) 

H  =%B  B  -  2J 

Fig.  2.  A  POSSIBLE  CONFIGURATION  OF  A  FINITE  SQUARE  LATTICE. 
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of  space  (clustering),  leaving  empty  space  in  the  rest  of 
the  container!]  11]  . 

For  ferromagnetism  the  total  internal  energy  is  the 
sum  of  the  interaction  energy  and  a  magnetic  energy  term 
-poBPi  for  each  lattice  point.  For  magnetic  systems,  Po  is 
a  characteristic  value  of  a  magnetic  moment  (such  as  the 
Bohr  magneton)  and  B  is  an  external  magnetic  field  (see 
figure  2).  For  other  physical  systems  this  term  -poBPi  may 
be  any  parameter  which  plays  the  role  of  a  "chemical 
potential"  in  determing  the  average  number  of  up  and  down 
spins,  or  average  composition  of  a  mixture,  or  average 
density  of  a  molecule-hole  system. 

For  ferromagnetic  materials  we  calculate  the  mean 
value  of  electron  magnetic  moment  as [10], 

M  =  NPotanhft 

where  N  =  number  of  atoms  per  unit  volume,  Po  is  the 
magnetic  moment  per  atom,  and  ft  represents  pBa/kBT  where  Ba 
is  the  mean  field  acting  on  the  atom,  and  kB T  is  the 
Boltzmann  energy  and  P  is  the  magnetic  moment  per  electron 
equal  to  q/2m  times  its  g-f actor,  times  its  angular  momentum 
J-> .  To  calculate  the  internal  energy  of  the  material,  we 
note  that  the  energy  of  an  electron  is  exactly  proportional 
to  the  magnetic  moment.  We  replace  Po  with  -PoB  to 
calculate  energy  and  B  can  be  written  as  (H  +  rM/'Soc2)  where 
H  is  the  magnetizing  field  and  r  has  been  called  the 
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"exchange"  force  and  is  due  to  the  exclusion  principle  in 
quantum  mechanics.  (Note:  we  are  using  r  for  the  "exchange 
force  in  this  discussion  to  avoid  confusion  with  J“ >  the 
angular  momentum.  In  the  rest  of  this  work  J  will  represent 
the  "exchange"  force).  Theoretical  predictions  for  the  value 
of  r  are  failures.  The  most  recent  calculations  of  the 
energy  between  the  two  electron  spins  in  iron-assuming  that 
the  interaction  is  a  direct  one  between  the  two  electrons  in 
neighboring  atoms-not  only  do  not  give  the  correct  value  but 
even  give  the  wrong  sign.  With  these  substitutions  the 
mean  energy  of  the  material  can  be  written  as, 

<U>  =  NPo  [H+rM/2£oc2  jtanhft 

The  "2"  is  inserted  to  correct  for  overcounting.  The  term 
rM/'Soc2  represents  interactions  of  all  possible  pairs  of 
atoms,  and  we  must  count  each  pair  only  once.  With  H  =  0  we 
can  rewrite  this  equation  as, 

M/Ms a t  =  tanh[(Tc/T)  (M/Maat  )] 
where  Meat  =  NjJ  and  Tc  =  PrMsat  /kSoc2  . 

The  solution  to  this  equation  is  shown  in  Fig. 3.  For 
the  energy  of  the  spontaneous  magnetization  below  the  Curie 
point,  we  can  set  H=0,  and  note  that  tanhft  =  M/Msat .  The 
mean  energy  is  proportional  to  M2 ,  and  can  be  written  as 
<0)av  —  -NMrM2  /  [  2*:  0  C2  Ms  a  t  ]  . 

If  we  plot  the  energy  due  to  the  magnetism  as  a  function  of 
temperature,  we  get  a  curve  which  is  the  negative  of  the 
square  of  curve  of  Fig. 3,  and  is  drawn  in  Fig. 4a.  If  we  were 
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34 


CRYSTALS  AS  A  FUNCTION  OF  TEMPERATURE  R.  FEYNMAN, 
THE  FEYNMAN  LECTURES  ON  PHYSICS.  2,  34-4,  (1964) 


u 


(a) 


Fig.  4.  THE  ENERGY  PER  UNIT  VOLUME  AND  SPECIFIC  HEAT  OF 
A  FERROMAGNETIC  CRYSTAL.  SOURCE:  SEE  Figure  3,  Page  34 


35 


■bo  measure  the  specific  heat  of  such  a  material  we  would 
obtain  a  curve  which  is  the  derivative  of  Fig. 4a,  and  is 
shown  in  Fig. 4b.  The  true  situation  is  more  complicated 
and  both  experiment  and  improved  theory  suggest  the  curve 
should  be  more  like  that  shown  in  Fig. 4c.  The  curve 
goes  higher  at  the  peak  and  falls  to  zero  somewhat  slowly. 

In  Fig. 5  several  specific  heat  curves  are  summarized 
for  a  two  dimensional  Lenz-Ising  model  calculated  with 
various  methods  of  approximation  [17].  Onsager’s  exact 
solution  for  B  =  0  is  shown  for  comparison[13] .  Approximate 
methods  [18,19]  have  attempted  to  reproduce  the  singularity 
in  the  specific  heat  and  other  thermodynamic  functions  at 
the  transition  point.  Methods  such  as  using  Pade’s 
approximants  [11]  to  extrapolate  series  expansions  of  the 
partition  functions  and  thermodynamic  properties  have  been 
used  in  the  past  twenty  years .  The  more  recent 
renormalization  group  work  has  given  the  most  accurate 
results  at  the  critical  point  [20]. We  are  using  a  different 
approach  to  understanding  the  behavior  of  the  model  over  a 
range  of  temperatures  and  magnetic  fields:  we  are 
expressing  the  free  energy  of  the  Lenz-Ising  model  using  the 
Morita  cluster  expansion  and  then  minimizing  the  free  energy 
expansion  to  obtain  the  equilibrium  macrostate  Po  and  then 
studying  the  resulting  approximation  to  the  thermal 
behavior . 
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Two-Dimensional 
sauare  lattice 


37 


Fig.  5.  SEVERAL  SPECIFIC  HEAT  CURVES  FOR  A  TWO  DIMENSIONAL  (SQUARE)  ISING  SPIN 
LATTICE  CALCULATED  WITH  VARIOUS  METHODS  OF  APPROXIMATION.  THE  EXACT 
TREATMENT  FOR  B=0  IS  ALSO  SHOWN.  R.  KUBO.  STATISTICAL  MECHANICS.  321,  (1965) 


Diagramatically  this  method  can  be  outlined  as  shown 
in  Fig. 6.  The  last  step  shown  in  Fig. 6,  to  "study  Fo " , 
means  to  study  these  behaviors : [21 , 22] 

"Study  Fo " 

U  =  F-TS  =  F+T(-5F/6T)b  i.e.  U  can  be 
obtained  from  F  and  its  derivatives. 


-M  =  ($F/<SB)t 

Cb  =  (JS/$T)b 

-S  =  (£F/<5'T)b 

x  =  ($M/$B)t 

where  M  is  the  magnetization 
S  is  the  entropy 
F  is  the  Helmholtz  free  energy 
T  is  the  temperature 
B  is  the  external  magnetic  field 
U  is  the  internal  energy 

Cb  is  the  specific  heat  at  constant  magnetic  field 
x  is  the  magnetization/particle 
In  this  work  we  also  use  3 ,  the  free  energy  in  reduced 
units : 

3  =  F/NJ, 

where  N  is  the  number  of  particles  in  the  system  and  J  is 
the  interaction  energy. 
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Fig.  6.  DIAGRAMATIC  OUTLINE  OF 
EXPANSION  OF  THE  FREE  ENERGY. 


CHAPTER  III 


EXACT  SOLUTION  OF  THE  LENZ-ISING 
MODEL  FOR  SMALL  SYSTEMS 

Exact  Solution  for  A  System  of  One  Particle 

The  Cluster  Expansion  Method  used  in  this  work 
involves  approximating  the  behavior  of  systems  of  many 
particles  as  a  concatination  of  simpler  systems:  systems  of 
one  particle,  systems  of  two  particles,  and  so  on.  In  order 
to  introduce  the  notations  and  methods  used  here  in  a  simple 
and  exact  context,  we  begin  by  studying  small  systems. 
Because  these  small  systems  can  be  studied  exactly,  their 
behavior  serves  as  a  measure  of  the  consequences  of  the 
approximations  used  in  the  Cluster  Expansion  Method. 

The  equilibrium  behavior  of  the  Lenz-Ising  system  of 
N  particles  for  fixed  temperature  T  and  magnetic  field  B  is 
given  by  Gibbs’  canonical  prescription: 

Z(T,B,  N)  =  EPexp(-flH(lO  ) 
is  the  canonical  partition  function, 

F(T, B, N)  =  -kBTlnZ 
is  the  equilibrium  free  energy, 

P<N)(p;T,B)  =  exp(-tfH(P))/Z, 
is  the  equlibrium  macrostate 

which  follows  exactly  from  minimization  of  F(P**  )  versis 

P(  N>  .  Since  there  are  2N  distinct  microstates,  Z(T,B,N)  has 
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2N  terms,  before  any  combining  of  similar  terms  is  effected; 
hence,  Z  can  be  explicitly  calculated  for  small  enough 
values  of  N.  When  N  =  10,  2N  =  1024;  direct  computation  of 
Z  is  feasible,  particularly  if  a  computer  is  used.  When  N  = 
20,  2N  ~  106 ;  direct  computation  of  Z  is  just  within  the 
bounds  of  feasibility  with  current  large  computers.  When 
N  =  30,  2N  ~  109;  direct  computation  of  Z  is  currently 
infeasible . 

In  order  to  introduce  the  notations  and  methods  used 
we  shall  discuss  the  cases  N  =  1 , 2 , 3 , and  4. 

Microstate : 

The  microstate  for  this  one  particle  system  is  a 
two -column 


at  any  instant,  the  system  is  either  spin-up  (P  =  +1)  or 
^pin-down  (P  =  -1). 

Energy: 

The  energy  H  is  H  =  H(P)  =  -poBM 


where , 


P  is  the  microstate, 

Po  is  the  magnetic  moment  of  the  particle, 
and  B  is  the  magnetic  field, 

The  exchange  term  -JPiPj  does  not  exist  here, 

since  there  can  be  no  pair  interaction  when  there  is  one 

particle.  The  system  is  an  ideal  paramagnet  which  is  a 

trivial  case  to  study. 

Macrostate : 

The  macrostate  is  the  probability  law  for  the 
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microstate , 


P  =  P(  l  )  (M  )  =  P(  1 ) 


where  "a"  is  the  probability  that  the  system  is  spin-up, and 

"b"  is  the  probability  that  the  system  is  spin-down.  Since 

"a"  and  ”b"  are  probabilities, 

a  £  0  constraint  1 

and  b  £  0 

and 


constraint  2 . 


a  +  b  =  1 


It  is  more  convenient  to  express  P(P)  in  terms  of  a  single 


number  controlled  by  a  single  constant,  mainly  the  average 


magnetization 


x  =  <P>  =  a  -  b  , 
-1  S  x  <  1  , 

with 

a  =  (1  +  x)/2 
b  =  (1  -  x)/2  ; 


with  -1  <  x  <  1 


Internal  Energy: 

The  internal  energy  of  this  spin  is 
U  =  <H>  =  Eip  =  -i  [H(P=+1)P<  l)  (P=+1)+H(P=-1)P<  i)  (P=-l)] 
=  H(P  =  +  l  )P<  1 )  (P  =1 )  +  H(P=-1)P(  i>  (P=-l) 

=  (-PoB)a  +  (PoB)b 
=  (-PoB)[*(l+x)]  +  (PoB)[^(l-x)] 

=  (PoB)[(-S)(l+x)  +  ^(1-x)] 

=  -(PoB)x  . 

Entropy: 

The  entropy  of  this  spin  is 
-S/kB  =  <lnP>  =  Eip  =  -i  P(  i )  (p  )lnP(  i )  (p  ) 
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=  P<  1)  (  - 1 )  lnP(  1)  (-1)  +  Pd)  ( 1 )  lnP(  l)  (1) 

=  a*ln(a)  +  b*ln(b) 

=  1+x)  ln^(  1+x)  +  1-x)  ln^(  1-x)  . 

Free  Energy: 

The  same  term  "free  energy"  is  used  for  two  distinct 
properties:  On  the  one  hand,  it  is  used  for  the  function 
F(P)  =  U(P)  -  TS(P) , 

where  P  is  any  macrostate  at  all,  not  necessarily  one 
corresponding  to  any  kind  of  equilibrium.  On  the  other 
hand,  it  is  used  for  the  function 
F(T,B)  =  F ( Pe  q ( T , B ) ) 

where  Peq  =  Peq(T,B)  is  the  unique  macrostate  describing  the 
system  in  equilibruim  at  a  temperature  T  and  in  a  magnetic 
field  B.  The  context  of  use  will  always  make  clear  which 
function  is  meant. 

The  free  energy  for  any  macrostate  P(x)  is 
F(x)  =  U(x)  -  TS(x) 

=  -PoBx  +  kB T[^( 1+x) ln^( 1+x) 

+  ^(l-x)ln^(l-x)]  . 

We  now  vary  x  to  locate  the  value  that  makes  F  a  minimum: 

0  =  (<5F/5-x) 

=  -Mo  B  +  kBT[^ln^(l+x)  -  ^ln^(l-x)]  . 

The  solution  of  this  equation  is 

xmin  =  tanh(M  o  B/kB  T)  =  tanh(ftPoB), 
where  d  -  1/kBT 
The  minimum  value  of  F  is 
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Fmin  =  F(xmin)  =  -P  o  BtanhfYJp  0  B] 

+  kBTO[l+tanh(|'ipoB)  ]  lnJ£[  l+tanh(ftP  o  B  )  ] 

+  l-tanhCtfPo  B)  ]  ln%[  l-tanh(ftPo  B)  ]  }  . 

Examination  of  (iS^F/Jx2)  shows  that  this  is  actually  a 

minimum  of  the  free  energy,  and  consequently  xmin  and  Fmin 

are  the  equilibrium  values: 

F  =  Fmi  n 

X  =  Xmi  n  . 

Other  Quantities: 

With  the  equilibrium  value  of  the  free  energy  F 
known,  we  can  find  the  equilibrium  value  of  other 
properties : 

S  =  -(JF/ffT)  =  -kB[*S(l±x)lnJS(l±x)]  . 

U  =  F  +  TS  =  <H>  =  -  PoBx  . 

Specific  Heat:  Cb  —  T(SS/ST!)  -  SU/6'T!  . 

Magnetization:  M  =  (SF/JB  =  NxPO  =  NPo  tanh(Po  B/kBT)  . 

Suspect ibility:  chi  =  (£M/<SB)t  . 

=  [NPB  2 /kB  T] sech2 (PB  B/kfi  T)  . 

For  PbB  <<  kBT, 
chi  =  Cp/T 
where  Cp  =  Np  b  2  /kB 
is  the  Curie  constant. 

A  System  of  Two  Particles 

We  derive  exactly  the  equilibrium  free  energy  F  and 
the  equilibrium  macrostate  P(  2 )  for  the  Lenz-Ising  model 
with  two  particles  as  a  function  of  the  external  magnetic 
field  B  and  the  temperature  T.  These  exact  results  will  be 
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obtained  using  the  approach  of  the  cluster  expansion  method 
in  order  to  illustrate  this  method  in  a  simple  case;  the 
method  of  Gibbs  ensemble  is  also  used.  Comparison  will  be 
made  between  these  exact  results  and  an  approximate  solution 
which  estimates  the  probability  distribution  as  the  product 
of  single  particle  probability  functions. 

Cluster  Expansion 

We  first  derive  the  exact  results  using  the  free 
energy  minimization  method.  The  free  energy  for  an 
arbitrary  macrostate  P< 2 )  is  given  by 
F(P<  2)  )  =  U(P<  2)  )  -  TS(P(  2)  ) 
where  U  is  the  internal  energy, 
and  S  is  the  entropy. 

For  a  specified  temperature  T  and  magnetic  field  B,  a  unique 
macrostate  P(2)min  minimizes  this  free  energy  F;  this 
macrostate  is  the  equilibrium  state  for  the  specified  T  and 
B,and  Fmi  n  =  F(P(2)min)  =  F(T,B)  contains  all  the  system’s 
equilibrium  behaviors. 

Microstate : 

The  microstate  of  the  system  is  M  =  (Ml  ,M2  )  where 
Mi  is  the  microstate  of  the  ith  particle. 

Energy : 

The  hamiltonian  H  is 

H  =  MoB(Mi+M2)  -  JM1M2, 

where 

B  is  the  external  magnetic  field 
Mo  is  the  magnetic  moment  of  each  particle 
and  J  is  the  exchange  integral  . 

Macrostate : 
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The  macrostate  of  the  system  is  P< 2 ) 


=  P(  2) 

1  1 

=  H 

l+xi  +x2  +y 

1  -1 

1+xi  -x2  -y 

-1  1 

1-xi  +X2  -y 

-1  -1 

1-xi  -x2  +y 

where  xi  =  <Pi  >  for  i  =  1,2 

and  y  =  <PiP2>  . 

Note  that  -1  i  xi  £  +1  and  that  -1  £  y  £  +1  . 

'  This  system  is  now  assumed  to  be  unchanged  under 
exchange  of  its  particles. 

Because  the  magnetic  field  B  is  uniform  in  space, 
so  that  it  acts  equally  on  each  particle,  and  because  all 
particles  have  the  same  properties,  then  all  equilibrium 
properties  must  be  unchanged  under  exchange  of  xl  and  x2 . 
Hence,  the  equlibrium  macrostate  can  depend  only  on  a  common 
value 


x  =  xi 
P(  2  )  =  \ 


-  X2  : 

l+2x+y 
1  -7 
1  "7 


l-2x+y 

and  this  simpler  expression  will  be  used  below. 
Internal  Energy: 

The  internal  energy  U  is 

U  =  <H>  =  -  (Po  B)  (Pi  +  P 2  )  -  JPi P 2 


=  -  2P  o  Bx  -  Jy. 


Entropy: 

The  entropy  S  is 


S  =  -kflEall  statesP<  2)  (p)lnP(  2)  (p) 

=  -kB[*(l+2x+y)ln(*)(.)  +2**( l-y)ln(S) ( . ) 
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+*(l-2x+y)ln(*)(. )] 


Free  Energy: 

The  free  energy  F  is, 

F  =  U  -  TS  =  - 2JJ  o  Bx  -  Jy  +  kB  T{  [34(  l+2x+y )  ]  In  [  .  ] 
+2[(*)<l-y)]ln[.]  +  [*(l-2x+y)]ln[.]}  . 

Note:  The  notation  [.]  means  the  quantity  which  precedes 

the  logarithm  is  repeated  as  the  argument  of  the  logarthim; 
that  is, 

[A] ln[ . ]  =  [A] ln[A]  . 

We  vary  x  and  y  in  the  free  energy  to  ( 1 )  locate 
those  values  of  x  and  y  which  minimize  F,  and  (2)  to  locate 
the  corresponding  minimum  value  of  F. 

(x-eqn. ) 

0  =  (SF/£x)y,  B,T 

=  -2P0B  +  kBT[(*)21n(*) (l+2x+y)  -  (%)21n(*) ( l-2x+y ) ] 
or  4rtPoB  =  ln[(l+2x+y)/(l-2x+y)]  . 

(y-eqn. ) 

0  =  (JF/$y)*,B,T 

=  -J  +  kBT[3$ln(%)  ( l  +  2x+y) 

-  2(*)ln(%)(l-y)  +  (S)ln(*)(l-2x+y)] 

( l+2x+y ) ( l-2x+y) 

or  4AJ  =  ln[  -  ]  . 

(l-y)<  2) 

These  are  two  simultaneous,  non-linear  equations  in  the  two 
unknowns  x  and  y.  Their  solution  gives  the  minimizing 
expressions  for  x  and  y: 

y  =  tanh{ft  J+(^)  ln[cosh(  2|3po  B)  ]  } 
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x  =  l+y)tanh(  2(1  Vo  B)  . 


(1) 


Note  that  this  solution  has  the  correct  behavior 
under  a  "mirror  reflection"  of  the  magnet.  Such  a 
reflection  causes  these  changes: 


x  --> 

7  --> 

B  --> 
T  --> 
J  — > 


-x 

7 

-B 

T 

J  ; 


These  are  compatible  with  equations  (1). 

When  J  =  0,  this  system  becomes  simply  a  pair  of  uncoupled 
particles,  that  is,  a  pair  of  the  systems  studied  in  the 
previous  section.  Indeed,  when  J  =  0,  equations  (1)  become 

y  =  x2  , 

i.e.,  the  particles  are  uncorrelated,  and 
x  =  tanhdPoB  , 

as  shown  in  the  previous  section. 

To  calculate  the  minimum  value  of  the  free  energy  we 
substitute  (x,y)min  into  F(x,y,;T,B),  obtaining  F(T,B). 

Gibbs  Ensemble: 

Following  the  suggestion  of  Dr.  J.  Goldman  we  next 
derive  the  same  exact  results  for  a  Lenz-Ising  model  of  two 
particles  using  the  Gibbs  ensemble  method. 

The  equilibrium  free  energy  is  given  by 


F  =  -kB  TlnZ 


where  the  partition  function  Z,is  expressed  by, 

Z  =  EPexp[-flH(P  )P  )  ]  ;  A  -  1/kB  T 

H(M)  is  the  energy  per  each  microstate 
=  -poB(Pi  +M2  )  -  JJJ  l  JJ  2  . 
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The  macro state  is 


P(M)  =  exp[-flH(P)]/Z  . 


For  the  Lenz-Ising  model  with  N  =  2,  the  microstates  are 


and  E(P ) 


Using  these  values,  we  calculate  Z  and  then  F: 

F  =  kBTln[4exp(r  )cosh(r+tf  J)  ] 
where  r  =  ^ln[cosh(  2/iJJo  B)  ]  . 


The  macrostate  is 


p<  /  +1  +1  \=  \ 

l+2x+y 

=  1/Z 

exp[-(3E(+l ,  +1 )  ] 

f  +1  -1  ) 

1  ~7 

exp[-flE(+l ,  -1 )  ] 

[-1+1 

1  -7 

exp[-(3E(  -1 ,  +1 )  ] 

\-l  -1  / 

l-2x+y 

exp[-GE(-l,-l)] 

*  / 

By  direct  comparison  we  find 


*(1  -  7)  =  exp[-UE(+l,-l)]/Z  ; 
hence,  y  =  tanh{(3  J+^ln[cosh(  2flpoB)  ]  }  , 

which  is  the  same  expression  for  y  as  obtained  previously 
using  the  cluster  expansion. 

Further,  by  direct  comparison,  we  find 

l+2x+y  =  4exp[i3  (  2Po  B+J )  ] /Z 
and  l-2x+y  =  4exp[fl  ( -2Po  B+J  )  ] /Z 

hence,  xg  =  {exp[(3  (  2Mo  B+J )  ] -exp[fl  ( -2Po  B+J )  ]  }/Z 
or, 


xg  =  xc l  uster  =  ^ ( 1+y ) tanh2ft P o B 
as  previously  obtained. 
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Magnetization : 


Magnetization  =  M  =  Po(Pi+P2)  =  2Pox 

=  (  2P  o  /2 )  ( l+y)  tanh2$Po  B 


Magnetic  Suspectibility : 

Susceptibility  =  chi  =  (6,M/i5'B)(  b  =  0 ) 

=  (  2P  0  /2 )  [  (i'  y/S  B )  (  B  =  0  )  ]tanh(2ftPoB) 

+(2Po /2)  (l+y)  (sech2  2|3po  B)  2i3Po 

=  (2flPo2  )  (1+tanhfJ  J) 

Comparison  of  Exact  P(  2 )  with  the  Macrostate 
for  Two  Uncorrelated  Paticles: 

This  system  has  two  particles,  and  the  behavior  of 
these  is  correlated  by  the  exchange  term  -JP1P2 .  To  measure 
the  extent  of  the  correlations  between  the  particles  one 
computes  the  macrostate  which 

(1)  is  consistent  with  the  exact  one-particle 


macrostates  P(  1 )  (  1 )  and  P<  1 )  (  2 )  , 

and  (2)  has  no  correlation  between  their  particles: 

[P(2)  (PlP2]uncorrelated  —  P(l)exactPlP(l)exactP2 

This  uncorrelated  macrostate  is 

[P(2)  (PlP2  )  ]unoor  r  el  at  ed  -  ^  1+X1  1+X2 

1-X1  1-X2 

(  1+X1  )  1+X2 

1-X2 


=  H 


1+X2 

(  1-X1  )  1-X2 


(  1+X1  )  (  1+X2  ) 
(  1+X1  )  (  1~X2  ) 


(  1-X1  )  (  1  +X2  ) 
(  1-X1  )  (  1-X2  ) 


1+X1  +X2  +X1  X2 
1+X1  “ X2  -XI  X2 


1  -XI  +X2  -XI  X2 
1-X1  —  X2  +X1  X2 


For  the  systems  of  interest  in  this  work,  xi  =  x2  =x;  hence, 
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l+2x+x2 

(1+X)2 

P(  2)  uncorrelated  =  \ 

1  -x2 

.  =  H  . 

(1+x) (1-x) 

1  -X2 

(l-x)(l+x) 

l-2x+x2 

.  ) 

(1-X)2 

k  J 

On  the  other  hand,  the  exact  case  is 
P(  2)  exact  -  H 


l+2x+y 
1  -7 
1  -7 
l-2x+y 


V  j 

If  the  two  particles  were  uncorrelated  (i.e.,  J<<kBT),  then 


P(  2  >  uncor  r  el  at  ed  =  P(  1  )  exact  *P(  O  exact  ,  and  SO 


y  -  <P  1  J->  2  >e  x  a  c  t  =  <P  1  >  <P  2  >  -  XI  X2  =  X2 
Hence  one  measure  of  the  extent  of  correlation  is  the  size 
of  the  difference  <PiP2>  -  <Pi  ><P2  =  y  -  x2 ;  this  difference 
is  zero  for  uncorrelated  particles. 

Other  Results: 

Other  results  for  uncorrelated  particles  are  as 

follows : 


Uuncor  r  el  at  ed  =  Hav  =  -  PoB(xi+X2)  -J<PlP2> 

=  PoB(xi+x2)  -Jyi  2 
=  -2PoBx  - Jx2  , 

-  (  S/kB  )uncorrelated  =  P<  2  )  (  ++  )  lnP(  2  )  (  ++  )  +P(  2  )  ( +- )  lnp<  2  )  ( +- ) 

+P<  2>  (-+)lnP<  2)  (-+)  +P<  2>  (  -  -  )  lnP<  2)  (--) 
=  Hi 1+x)2 ln( . )  +  2(H) (1-x2 )ln^(l-x2 ) 

+Hi 1-x)2 ln( .  ) 

=  2(^)(l+x)2ln^(l+x)+2(?4)(l-x2  )ln^(l-x2  ) 

+2(*)(l-x)21n*S(l-x) 

=  Hi 1+x)2  +  Hi l~x2 ) ln^( 1+x) 

+^(l-x)2  +  J$(l-x2  )ln^(l-x) 
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=  2[JS(l+x)ln)*(l+x)  +  ^(l-x)ln^(l-x)  )  . 

=  -(Si  +  S2  )/kB  , 

or  S(  2  )  uncor  t  el  ated  =  Sl(^)  +  S2  (  1  ) 

The  approximate  free  energy  is  given  by, 

Fun  correlated  — - 2M  o  Bx— J x2  +2kT [ ^ ( 1 +X ) ln^ ( 1 +x ) ( 1 —  X ) ln^( 1  —  X ) ] 
=  F(  1 )  l  +  F(  1 )  2  -  Jx2  .  (  2 ) 

This  equation  for  the  free  energy  can  be  minimized  versis 
changes  in  x: 

(x-eqn)  0  =  £F/£x  =  -2jJoB  -2Jx  +2kT [^ln^(  1+x)  -^ln^(l-x)] 
(JJoB+Jx)kT  =  JSln[(l+x)/(l+x)]  =  tanh(-l)x  , 
or  x  =  tanh[(PoB)/kBT  +  (J/kBT)x] 

=  tanh((ipoB  +$Jx)  (3) 

where  (1  =  1/kBT. 

When  J  =  0,  the  two  particles  are  independent,  and  hence 
they  are  uncorrelated;  intensive  properties  are  exactly  the 
sum  as  in  the  N=1  case  studied  above,  while  extensive 
properties  are  exactly  the  sum  of  the  extensive  properties 
of  the  individual  particles.  When  J  >  0,  then  the  two 
particles  interact,  and  there  is  correlation.  However, 
one  can  retain  the  interaction  but  ignore  the  correlation, 

by,  for  example,  using  Eq.  (2)  for  the  free  energy. 

The  result  for  x,  Eq.  (3),  is  an  approximation  to  the 
correct  results;  this  approximate  result  for  x  can  be  used 
to  obtain  approximate  results  for  other  quantities, 
including  U,  S,  and  F. 
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Summary 

The  same  exact,  results  have  been  derived  for  N  =  2 
Lenz-Ising  model  using  either  the  cluster  expansion  or  the 
Gibbs  ensemble  method.  The  Gibbs  ensemble  can  only  be  used 
for  exact  calculations,  whereas  the  cluster  expansion  can  be 
used  for  either  exact  or  approximate  calculations.  We 
found  that  for  the  minimum  free  energy  the  agreement 
between  exact  and  approximate  results  was  good. 

It  is  noted  that  for  N  =  2  the  results  for  B  =  0  are 
quite  different  than  those  for  B  =/  0,  whether  the  results 
are  obtained  by  exact  or  approximate  methods.  Since  many  of 
the  published  calculations  for  the  Lenz-Ising  model  are  for 
the  case  B  =  0  then  it  will  be  important  to  see  if  this 
tendency  for  non-zero  B-field  results  to  be  different  than 
zero  B-field  results  continues  for  larger  number  of 
particles . 

A  System  of  Three  Particles 
We  derive  the  exact  results, for  a  Lenz-Ising  model 
with  only  three  particles  using  the  cluster  expansion. 

There  are  several  possible  configurations  for  three 
particles;  Case  I-  three  particles  on  a  straight  line,  Case 
II-  an  equilateral  triangle  (all  J’s  equal),  Case  Ill-a 
right  triangle.  We  begin  with  Case  II,  and  consider  an 
equilateral  triangle. 

Microstate : 

The  microstate  of  this  system  is  P  -  (Pi P2P 3  )  where 
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Pi  is  the  microstate  of  the  ith  particle. 

Energy: 

The  hamiltonian  for  this  system  is, 

H  =  ~PoB(Pl  +P2 +P3  )  -  J(PlP2+P2P3+PlP3  ) 

where 

Pi  ,  P 2 ,  P 3  are  the  1-particle  microstates  for 
particle  1,2,3  respectively, 

B  is  the  external  magnetic  field 
and  J  is  the  interaction  term  for  pairs  of 
particles . 

Macrostate : 

There  are  three  particles  and  initially  we  do  not 
treat  them  as  similar.  Later  we  will  let  them  be  the  same. 
The  microstates  are: 


Pi 

1 

1 

1 

1 

-1 

-1 

-1 

-1 


P  2 
1 
1 

-1 

-1 

1 

1 

-1 

-1 


P  3 
1 

-1 

1 

-1 

1 

-1 

1 

-1 


This  is  an  exact 
configuration . 


Each  macrostate  is  a  probability  law  for  these  8 
microstates.  The  macrostate  P(  3 )  is 

P(  3)  (P1P2P3  )  =  P(  3) 


Pi 

P  2 

P  3 

3)  1 

1 

1 

s 

1 

1 

-1 

h 

1 

-1 

1 

i 

1 

-1 

-1  = 

3 

-1 

1 

1 

k 

-1 

1 

-1 

1 

-1 

-1 

1 

m 

-1 

-1 

-1 

n 

constraints : 

g,h, 

i  ?  3  t  k  )  1 : 

,m,n 
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The  quantities  g  ->  n,  are  just  numbers.  We  can  either 
describe  probability  functions  in  terms  of  numbers  such  as 
these,  or  in  terms  of  physical  quantities  x,y,z.  The 
advantage  of  working  with  g  ->  n  is  that  they  do  not  mean 
anything.  They  are  just  numbers,  and  only  mean  the 
probability  of  being  up  or  the  probability  of  being  down. 
For  example,  the  probability  for  the  microstate  with  all 


three  particles 

being  up, 

is 

given 

by  "g 

it  i 

The  physical 

quantities 

,  such  as 

"x"  , 

mean 

something 

-  it 

is  the 

probability  of 

spin 

up  minus 

the  probability 

of  spin  down 

i.e.  x  =  p(~)  - 

P(  ) 

• 

The 

macrostate 

for 

this 

three  particle 

system  is, 

P<  3)  (PI  ,P2 

,P3  ) 

f 

XI 

X2 

X3 

yi  2 

yi  3 

y2  3 

Z1  2  3'' 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

-1 

1 

-1 

-1 

-1 

1 

1 

-1 

1 

-1 

1 

-1 

-1 

=  1/8  < 

1 

1 

-1 

-1 

-1 

-1 

1 

1 

► 

1 

-1 

1 

1 

-1 

-1 

1 

-1 

1 

-1 

1 

-1 

-1 

1 

-1 

1 

1 

-1 

-1 

1 

1 

-1 

-1 

1 

,  1 

-1 

-1 

-1 

1 

1 

1 

"I  , 

where  x,y,z,  have  been  previously  defined.  Note  that  the 
individual  y  entries  are  products  of  x’s  (yi2  is  the  product 
of  xi  ,  X2  etc.  )  and  zi2  3  is  the  product  of  xi  ,  X2  ,  X3  .  For 
example  in  the  second  row  y2  3  =  (x2)(x3)  =  ( 1 )  ( —  1 )  =  -1,  and 
in  the  third  row  zi23  =  (1)(-1)(1)  =  -1. 

How  can  we  check  this  arrangement?  Use  normalization 
and  add  up  the  g,h,...m,n’s.  For  example, 


g  =  1/8 


x, s  y’s  z ’ s 

0  0  0 


1/8(8) 
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For  the  Lenz-Ising  model  these  three  particles  will  be 


similar,  and  we  can  add  up  the  x’s  and  the  y’s.  For 
example,  3x  =  xi  +  X2  +  X3  and  all  three  spins  are  pointing 
up.  Then  P< 3 )  can  be  written  as, 

/ 


1 

+  3x 

+  3y 

+  Z 

1 

+x 

-y 

-z 

1 

+x 

“7 

“Z 

P<  3) 

=  1/8 

< 

1 

-x 

“7 

+z 

1 

+x 

“7 

-z 

1 

-x 

"7 

+z 

1 

-x 

“7 

+  2 

1 

-3x 

+  37 

-z 

L 

> 

Remember 

that  the 

X  < 

and  7 

entries 

above  are 

3y  =  yi  2 

+  yi  3  + 

y2  3 

(the 

sum  for  all  three 

value  for  z  is  a  product  of  the  three  x  values. 

Internal  Energy: 

The  internal  energy  is  the  the  average  value  of  the 
hamiltonian,  i.e.  the  expectation  value  of  H, 


U  =  <H>  =  -  PoB[<Pi  >+<P2  >  +  <P3  >] 

-  J[<PlP2 >+<P2P3 >+<PlP3 >] 

Recall  the  N=1  case.  We  defined  the  average  value  or  the 

expectation  of  the  spin-state,  P ,  as 

x  =  <P  >  =  Ep-i  1  P<  i  >  (P  ) 

=  (p=+l)  (P<  l)  (+1)  )  +  (P=-1)(P<  i)  (-1)  ) 

We  recognize  symmetries  in  the  crystal  lattice, 

<p  1  >  =  <p  2  >  =  <P  3  >  =  X 

the  average  value  of  P ,  i.e.  the  average  value  of  the  spin 
state ,  and 

<P  1  P  2  >  =  <P  2  P  3  >  =  <P  1  P  3  >  =  y 
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the  average  value  of  the  product  of  P . 

The  internal  energy  can  then  be  written  as, 

U  =  -3PoBx  -  3Jy 

where  x  is  the  average  value  of  the  magnetization  of  the 
particle.  It  is  the  probability  of  spin  up  minus  the 
probability  of  spin  down:  i.e.  x=p(~)-p(  ).  The  quantity  y 

is  the  (average  value)  probability  of  finding  a  pair  being 
parallel . 

Entropy: 

The  entropy  for  this  system  is, 

-S/kB  =  <lnPO)> 

-S/kB  =  [E  ,all  states]P(  3>  (P1P2P3  )  lnP(  3)  (P1P2P3  )  . 

The  exact  configuration  for  the  microstatees  is  a  3x8  array 
(see  p.54).  We  wrote  the  macrostate  in  terms  of  the 
numbers  g  -->  n.  To  write  this  macrostate  in  the  form  we 
use  with  xi,x2,yi2,z  we  need  to  evaluate  the  numbers  g  -->  n 
as  we  evaluated  the  analogous  numbers  c,d,e,f  for  P(  2 )  .  For 
the  N=3  case  the  new  quantity  is  z  =  <PiP2P3>  i.e.  the 
product  of  all  three  microstates.  The  quantities  x  and  y 
are  the  same  as  in  macrostate  P<  2  >  . 

REDUCING 

Is  it  the  case  that  if  I  sum  P(  3)  over  any  one  of 
these  three  particles  I  get  the  P(  2 )  of  the  other  2 
particles.  We  do  this  in  the  following  way,  by  extracting 
the  behavior  of  2  particles  from  the  known  behavior  of  3 
particles.  We  do  this  by  adding  up  over  all  possible  states 
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of  the  third  particle.  This  is  known  as  reducing,  (reducing 
the  marginal  probability) 


P(  2  )  (JJ1P2)  =  EP3=-1  1  P<  3  )  (JJHJ2M3  ) 
P(  2)  ( JJ  1 JJ  3  )  =  S)J2  =  -1  1  P(  3)  (Pi  P2P  3  ) 
P(2)(p2P3)  =  IPl  =  -ll  ]P<  3)  (P1P2P3  ) 


P(  2)  (pip  2  ) 

11  111  11-1 

=  P(2)  1-1  =  P(3)  -1  11  +  P(3)  1  -1  -1 

-11  -111  -11-1 

-1  -1  -1  -1  1  -1  -1  -1 

P3=l  P  3  = - 1 

P(  2)  (p  1  p  2  ) 

(1  ID  (11  -1) 

=  1/8 ( l+3x+3y+z )  +1/8 ( 1+x-y-z )  — > 1/8 ( 2+4x+2y ) 

(-1  1  1)  (1  -1  -1) 

1/8 ( 1+x-y-z )  +1/8 ( 1-x-y+z )  -->  1/8(2  -2y) 

(-1  1  1)  (-1  1  -1) 

l/8(l+x-y-z)  +1/8 ( 1-x-y+z )  -->  1/8(2  -2y) 

(-1-1  1)  (-1-1-1) 


l/8(l-x-y+z)  +1/8 ( l-3x+3y-z )  -->  1/8(2  -4x+2y) 

=  \  l+2x+y 
1  "  7 
1  -  7 
1-2+y 

This  is  the  same  P< 2 )  as  derived  earlier  -  this  time  by  a 
reduction  method. 

Entropy : 

-  S/k  =  SP1P2P3P<  3)  lnP<  3) 

=  ( 1/8 ) ( l+3x+3y+z ) ln( 1/8 ) ( l+3x+3y+z ) 

+(3/8) (1+x-y-z )ln( 1/8) (1+x-y-z) 

+ ( 3/8 ) ( 1-x-y+z ) ln( 1/8 ) ( 1-x-y+z ) 

+ ( 1/8 ) ( l-3x+3y-z ) In (1/8) (l-3x+3y-z) 

Energy : 

F  =  U  -  TS 

=  -3P  0  Bx  -  3 Jy 

+  (kB  T/8 ) [ ( l+3x+3y+z ) ln( 1/8 ) ( l+3x+3y+z ) 

+3{ ( 1+x-y-z ) ln( 1/8 ) ( 1+x-y-z ) } 

+3 ( 1-x-y+z ) In ( 1/8 ) ( i-x-y+z ) 

+ ( l-3x+3y-z ) In (1/8) (l-3x+3y-z)] 

We  now  MINIMIZE 

0  =  6'  F  /6'  x  =  -3P  0  B  +  (kBT/8)  [3  In  (1/8)  (l  +  3x+3y+z) 

+31n( 1/8 ) ( 1+x-y-z ) 

-31n( 1/8 ) (1-x-y+z) 

-31n( 1/8 ) ( l-3x+3y-z) ] 
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0  =  SF/Sy 

-  -3J  +  (kBT/8) [3 In (1/8) ( l+3x+3y+z ) 
-31n(l/8) (1+x-y-z) 

-31n( 1/8 ) (1-x-y+z ) 

+3 In ( 1/8 ) ( l-3x+3y-z ) ] 

0  =  SF/Sz 

=  (kBT)[ln(l/8)(l+3x+3y+z) 

-31n(l/8) (1+x-y-z) 

31n( 1/8 ) ( 1-x-y+z ) 

-ln( 1/8 ) ( l-3x+3y-z ) ] 

Suppose  B=0 ,  so  we  can  study  (x,y,z)  with  B=0 


[x  eqn] 

( l+3x+3y+z ) ( 1+x-y-z ) 

-  =  1 

( l-3x+3y-z) ( 1-x-y+z) 

— >  (l+3x+3y+z) (1+x-y-z)  =  ( l-3x+3y-z ) ( 1-x-y+z ) 


[z  eqn] 

( l+3x+3y+z ) ( 1-x-y+z ) 3 

- -  1 

( l-3x+3y-z ) ( 1+x-y-z ) 3 

— >  ( l+3x+3y+z )( 1-x-y+z )3  =  ( l-3x+3y-z )( 1+x-y-z ) 
From  [x] ,  we  know 


From  [z] , 

(l+3x+3y+z)  (1-x-y+z) 

(l-3x+3y-z)  (1+x-y-z) 

we  know 

(l+3x+3y+z)  (1+x-y-z)3 

(l-3x+3y-z)  (1-x-y+z)3 

L  Hence 

(1-x-y+z)  (1+x-y-z)3 

(1+x-y-z)  (1-x-y+z)3 

(1-x-y+z)4  =  (1+x-y-z)4 

1-x-y+z  =  1+x-y-z 
-x+z  =  x-z 
[x  =  z] 

Return  to  [x] ,  with  (x  =  z) 

(l+4x+3y)/(l-4x+3y)  =  (l-y)/(l-y)  =  1 
-->  l+4x+3y  =  l-4x+3y 

— >  [x  =  0] 
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Hence,  (B  =  0)  -->  [x  =  z  =  0] 

We  continue,  B  =  0  with  [y  eqn] 

3J  =  2kBT/8[31n(l/8) ( l+3y ) -31n( 1/8 ) (1-y)] 

—  >  J  =  (ka  T/4 ) ln[ ( l+3y ) /( 1-y ) ] 

4J/kB  T  =  ln[  ( l  +  3y )/( 1-y )  ]  ot  =  4J/kBT  =  4rtJ 

*  =  ln[ ( l+3y)/( 1-y ) ]  — >  exp(a)  =  (l+3y)/(l-y) 


-->  exp(a)  -  exp(a)y  =  l  +  3y 

-->  exp(<x)  -  1  =  y[3  +  exp(a)] 

y  =  [exp(a)  -  l]/[3+exp(oc )  ] 

y  =  [exp(4flJ)  -  l]/[exp(4ftJ)  +  3] 


Case  I 

We  have  been  examing  case  II  for  triangle 
arrangements.  There  are  two  other  possibilities  for  a  three 
particle  Lenz-Ising  lattice  models:  the  triangle  and 
straight  line.  We  now  examine  case  I  for  3  particles  in  a 
straight  line. 

For  a  straight  line  we  have  Ji2=J23=J  and  Ji 3  ~  0 
Then  yi2=y2  3=yi  and  yi  3 =y2 .  The  microstates  are  as  before, 


P(  3) 


111 
11-1 
1-11 
1  -1  -1 
-1  1  1  [ 
-1  1  -1 
-1  -1  1 
-1  -1  -1 


The  new  macrostate 

/ 


P<  3)  =  1/8 


y 

is, 


1 

+ 

3x 

+ 

2yi 

+ 

y2 

1 

+ 

X 

- 

y2 

1 

+ 

X 

- 

2yi 

+ 

y2 

1 

- 

X 

- 

y2 

1 

+ 

X 

- 

y2 

1 

- 

X 

- 

2yi 

+ 

y2 

1 

- 

X 

- 

y2 

1 

- 

3x 

+ 

2yi 

+ 

y2 

z 

z 

z 

z 

z 

z 

z 

z 
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Free  Energy: 

U  =  -3(->oBx  -  2Jyi 
F  =  -3MoBx  -  2Jyi 


+kB  T{  ( 1/8  )  ( 1 

+ 

3x 

+  2yi 

+ 

72 

+  z ) ln( 1/8 ) ( . ) 

+ ( 2/8 ) ( 1 

+ 

X 

-  y2 

-  z ) ln( 1/8 ) ( . ) 

+(1/8)(1 

+ 

X 

-  2yi 

72 

-  z ) ln( 1/8 ) ( . ) 

+( 2/8 ) ( 1 

- 

X 

-  y2 

+  z ) ln( 1/8 ) ( . ) 

+(1/8)(1 

- 

X 

-  2yi 

+ 

72 

+  z) In (1/8) ( . ) 

+(1/8) (1 

— 

3x 

+  2yi 

+ 

72 

-  z ) ln( 1/8) ( . ) 

1:  Suppose  B 

0  . 

Then 

X 

= 

z  =  0 . 

0  =  SF/Syi  =  -2J  +  (kB  T/8  )  [ 21n(  1/8  )  ( l  +  2yi  +y2  ) 
[yi]  -2 In (1/8)  (l-2yi  +y2  ) 

-21n( 1/8 ) ( l-2yi +y2 ) 
+21n( 1/8 ) ( l+2yi +y2 )] 

0  -  SF/Sy2  -  (kT/8)[ln(l/8)(l  +  2yi+y2) 

-21n( 1/8 ) ( l-y2 ) 

+  ln( 1/8 ) ( l-2yi +y2 ) 

[y2  ]  -21n( 1/8 ) ( l-y2 ) 

+  ln( 1/8 ) ( l-2yi +y2 ) 

+  ln( 1/8 ) ( l+2yi +y2 )] 

Solve  [yi ]  and  [y2 ] : 

[yi  ]  J  =  (kB  T/4 )  ln[  ( l+2yi +y2  )  /  ( l-2yi +y2  )  ] 

—  >  4J/kBT  =  ln[  ( l+2yi +y2  )/( l-2yi +y2  )  ] 

[72  ] 

0  =  ln[  ( l+2yi +y2  )2  ( l-2yi +y2  )2 /( l-y2  )4  ] 
-->  1  =  ( l  +  2yi +y2  )  ( l-2yi +y2  )/( l-y2  )2 

( l+y2  +2yi  )  ( l+y2 -2yi  )  =  (l-y2)2 
(l+y2)2  -  4yi 2  =  (l-y2)2 
1  +  2y2  +  y22  -4yi 2  =  1  -  2y2  +  y22 

4y2  =  4yi 2 
[y2  =  yi  2  ] 


Substitute  this  result  into  [yi  ] 

[yi  ]  4J/kBT  =  ln[  ( l  +  2yi +yi  2  )  /  ( l-2yi +yi  2  )  ] 
2(J/kBT)  =  ln[  ( 1+yi  )/(l-yi  )] 

Results : 
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yi  =  tanh(ftJ) 

72  =  [tanh(l3J)]2 
For  B  =  0  and 
x  =  0 
z  =  0 

straight  line  tanh(x)  =  [l-exp( -2x) ]/[l+exp(-2x) ] 

case 


A  System  of  Four  Particles 

This  will  be  a  4-bod7  exact  calculation  for  an  equi- 
ever7thing  p7ramid.  In  these  exact  calculations  it  will  be 
the  first  time  that  we  have  encountered  a  three  dimensional 
Lenz-Ising  model.  For  this  P7ramid  all  sides,  faces  and 
angles  are  equal . 

We  have , 

x  =  <Pl  >  =  <P  2  >  =  <P  3  >  =  <P  4  > 

Which  is  the  average  behavior  of  each  particle  separate^. 

7  =  <P  1  P  2  >  =  <P  1 P  3  >  =  <P  1  P  4  > 

=  <P  2  P  3  >  =  <P  2  P  4  >  =  <P  3  P  4  > 

i.e.  the  average  correlation  for  each  pair  of  particles. 

z  =  <P  1 P  2  P  3  >  =  <P  1 P  2  P  4  >  =  <P  1 P  3  P  4  > 

=  <P  2  P  3  P  4  > 


W  =  <PlP2P3P4> 

Macro state : 

The  behavior  of  a  single  particle  is  given  b7  the 
macrostate  P<  1 ) 


1+x 

=  p<  1 ) 

1 

1-x 

-1 

The  behavior  of  one  particle  with  another  is  given  b7 


the  macrostate  P<  2 ) 

“  “ 

1+2x+7 

1  1 

P(  2)  (p  p  >  )  =  \ 

1  "  7 

1  -  7 
1-2x+7 

=  P(  2) 

1  -1 
-1  1 
-1  -1 

. 
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P<  2)  (p  p  ’  )  =  P(  2)  (Pi  !  Pi  2  )  . 


This  is  an  alternative 

way  of  writing  P(  2  >  by  using 

the  Morita  notation. 


For  three  particles,  the  macrostate  P< 3 )  is 


P<  3)  (M  ,P  '  ,M  ’  ’  )  =  1/8 


l+3x+3y+z 

1+x-y-z 

1+x-y-z 

1-x-y+z 

1+x-y-z 

1-x-y+z 

1-x-y+z 

l-3x+3y-z 


=  P<  3  ) 


111 
1  1-1 
1-1  1 
1-1-1 
-111 
-1  1-1 
-1-1  1 
-1-1-1 
microstates 


For  four  particles,  the  macrostate  P<4)  is 


=  P(  4) 


»  M  >  > 

)  “ 

,P  ’  ’ 

’  ) 

1 

1 

1 

1 

■p  l+4x+6y+4z+w  _ 

1 

1 

1 

-1 

l+2x+0y-2z-w 

1 

1 

-1 

1 

l+2x+0y-2z-w 

1 

1 

-1 

-1 

1  l+0x-2y+0z+w 

1 

-1 

1 

1 

l+2x+0y-2z-w 

1 

-1 

1 

-1 

l+0x-2y+0z+w 

1 

-1 

-1 

1 

l+0x-2y+0z+w 

1 

-1 

-1 

-1 

=  1/16 

l-2x+0y+2z-w 

'  "I 

1 

1 

1 

l+2x+0y-2z-w 

-1 

1 

1 

-1 

l+0x-2y+0z+w 

"I 

1 

-1 

1 

l+0x-2y+0z+w 

;  -1 

1 

-1 

-1 

l-2x+0y+2z-w 

-1 

-1 

1 

1 

l+0x-2y+0z+w 

;  -1 

-1 

1 

-1 

l-2x+0y+2z-w 

-1 

-1 

-1 

1 

l-2x+0y+2z-w 

-1 

-1 

-1 

-1  J 

—  l-4x+6y-4z+w  “ 

microstates  of  the 
four  particles. 


16  0  0 


Internal  Energy: 

U  =  <H>  =  - (Po  B) ( 4x)  -  ( J ) ( 6 ) ( y ) 

U  =  -  4PoBx  -  6Jy 

Entropy: 

-  S/kB  =  IPP<  4>  lnP<  4> 

=  (1/16) (l+4x+6y+4z+w)ln(l/16) ( .  ) 
+4(1/16) (l  +  2x-2z-w)ln( 1/16) ( .  ) 
+6(1/16) (l-2y+w)ln( 1/16) ( .  ) 
+4(1/16) (l-2x+2z-w)ln( . ) 

+  ( 1/16 )(l-4x+6y-4z+w)ln( 1/16 )( .  ) 
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Free  Energy: 


F  =  -  4MoBx  -  6Jy 

+  (kBT/16)[( l+4x+6y+4z+w) In ( 1/16 ) ( . ) 

+4 ( l+2x-2z-w) ln( 1/16 ) ( . ) 

+6 ( l-2y+w) In ( 1/16 ) ( . ) 

+4 ( l-2x+2z-w) ln( 1/16 ) ( . ) 

+( l-4x+6y-4z+w) ln( 1/16 ) ( . ) ] 


Minimization: 

[x]  0  =  SF/Sx 

=  -4PoB  +  (kBT/16) [4 In (1/16) ( l+4x+6y+4z+w) 

+8 In (1/16) ( l+2x-2z-w) 

-8 In (1/16) ( l-2x+2z-w 
-41n( 1/16 ) ( l-4x+6y-4z+w) ] 

( l+4x+6y+4z+w) (l+2x-2z-w)2 

—  >  4)J o B/kB T  =  %ln[  -  ] 

( l-4x+6y-4z+w) (l-2x+2z-w)2 

[y]  0  =  SF/Sy 

-  -6J  +  (kBT/16) [61n(l/16) (l+4x+6y+4z+w) 

-12 In (1/16) ( l-2y+w) 

+61n( 1/16 ) ( 1 -4x+6y-4z+w ) ] 

( l+4x+6y+4z+w) ( l-4x+6y-4z+w) 

—  >  16J/kB T  =  ln[  -  ] 

(l-2y+w)2 

[z]  0  =  SF/Sz 

=  (kBT/16) [41n(l/16) (l+4x+6y+z+w) 

-8 In (1/16) ( l+2x-2z-w) 

+81n(l/16) (l-2x+2z-w) 

-4 In ( 1/16 ) ( l-4x+6y-4z+w) ] 

( l+4x+6y+4z+w) ( l-2x+2z-w)2 

—  >  0  =  ln[  - ] 

( l-4x+6y-4z+w) ( l+2x-2z-w)2 

( l+4x+6y+4z+w) ( l-2x+2z-w)2 

—  >  1  =  - 

( l-4x+6y-4z+w) ( l+2x-2z-w)2 

[w]  0  =  SF/Sv 

-  (kBT/16) [ln(l/16) (l+4x+6y+4z+w) 

-4 In (1/16) (l+2x-2z-w) 

+61n( 1/16 ) ( l-2y+w) 

-41n(l/16) (l-2x+2z-w) 

+ln( 1/16 ) ( l-4x+6y-4z+w) ] 

( l+4x+6y+4z+w) ( l-4x+6y-4z+w) ( l-2y+w)6 

-->  0  =  ln[  -  ] 

(l+2x-2z-w)4 (l-2x+2z-w)4 
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( l+4x+6y+4z+w) ( l-4x+6y-4z+w) ( l-2y+w)6 

—  >  i  -  - 

( l+2x-2z-w) 4 (l-2x+2z-w)4 


Suppose  B  =  0 . 

[x]  ( l+4x+6y+4z+w) ( l+2x-2z-w)2 

- -  ! 

( l-4x+6y-4z+w) ( l-2x+2z-w)2 

[z]  ( l+4x+6y+4z+w) ( l-2x+2z-w)2 

- =  ! 

( l-4x+6y-4z+w) ( l+2x-2z-w)2 

( l+4x+6y+4z+w)  (l-2x+2z-w)2 

- - -  ( from  [x]  ) 

( l-4x+6y-4z+w)  (l+2x-2z-w)2 


( l+2x-2z-w)2 


ditto 

( l-2x+2z-w)2 

(from  [z] ) 

—  > 

(l+2x-2z-w)  = 

( l-2x+2z-w) 

—  > 

4x  = 

4z 

--> 

C  x  = 

z  ] 

Using  this  result  in  the  [  x  ]  equation: 


l+8x+6y+w 

- =  1 

l-8x+6y+w 

-->  l+8x+6y+w  =  l-8x+6y+w 

-->  16x  =  0 

—  >  x  =  0 

Hence,  (  B  =  0  )  -->  (  x  =  z  =  0  ) 

Using  x  =  z  =  0,  return  to  the  original  equation  for  [x] . 

( l+6y+w) ( 1-w)2 

[x]  16Po  B/kB  T  =  In  [ - ] 

( l+6y+w) (1-w)2 

0  =  lnl  since  B=0 

0=0 
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i.e.,  the  x-eqn  is  automaticaly  satisfied 


[7l 

( l  +  6y+w)2 

16ft  J  =  In  [ - ] 

( l-2y+w)2 

l+6y+w 

8GJ  =  In  ( - ) 

l-2y+w 

(4) 

[z] 

is  automatically  satisfied,  (  1=1  ) 

[w] 

( l+6y+w)2 ( l-2y+w)6  =  (1-w)® 

—  > 

( l+6y+w) ( l-2y+w)3  =  (l-w)4 

— > 

(  1  —  W  )  4 

l+6y+w  =  -  - 

(l-2y+w)3 

(5) 

The  above  two  equations  (4)  and  (5)  for  y  and  w 
are  two  simultaneous  equations  in  two  unknowns. 
Substitute  for  l+6y+w  in  the  eqn.  for  [  y  ]. 

(1-w)4 

8(iJ  =  In  [ - ] 

(l-2y+w)4 
1-w ' 

2ft  J  =  In  ( - ) 

l-2y 

-->  o<  =  exp(2i3)  =  ( l-w)/( l-2y+w) 

—  >  a  -  2<xy  +  otw  =  1  -  w 

—  >  a  -  1  -  2oty  +  (a  +  l)w  =  0 

—  >  (<x  +  l)w=l-a  +  2ay 


-->  w  =  (1  -  a )/( 1  +  a)  +  (2ay)/(l  -  a) 


Hence,  the  Free  Energy  is  for  B  =  0 

F  =  -6Jy  +  (kB  T/16 ) [2 ( l+6y+w) In ( 1/16 ) ( . ) 

+  8 ( 1-w) ln( 1/16 ) (  ) 

+  6 ( l-2y+w) In ( 1/16 ) ( . )] 
with  w  as  above . 
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Write , 


a  =  |3F  =  -  6(3  Jy  +  ( 1/8  )  [  ( l+6y+w) ln(  1/16 )(.  ) 

+4 ( 1-w) ln( 1/16 ) ( . ) 

+3 ( l-2y+w) ln( 1/16 ) ( .  )] 

In  the  Morita  1-cluster  approximation, 

y  =  <  PP’>  »  <  P  ><  P’>  =  <  P  X  P  >  =  x2  =  0;  B  =  0 

z  =  <  p  p  ’  p  ’  ’  >  s  <jj  >  <  p  ’  >  <  p  ‘  *  >  =  X3 

w  =  <  p  p*p»»p*»»>  z  <p><p  ’  ><p  ’  '  ><p  *  *  '  >  =  X4 

In  the  Morita  2-cluster  approximation, 

y  =  <  PP ’ >  is  treated  exactly. 

z  =  <  pp  »p  >  ’  >  =  <  p  ><p  >p  >  >  >  =  Xy  =  0 ;  B  =  0 

w  =  <pp’p,,p’’,>  =  <pp  ’  ><p  ’ ’p  ’  ’  ’  >  =  y2  . 

In  the  Morita  3-cluster  approximation, 
y,z  are  treated  exactly, 

w  =  <  PP’X  p  ’  ’p  ’  ’  ’  >  +  <  p  ><  p  p ’p  ’  ’  > 
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CHAPTER  IV 


APPLICATION  OF  THE  CLUSTER  VARIATIONAL 
METHOD  TO  THE  LENZ-ISING  MODEL 

First  Order  Approximation:  Mean  Field  Theory 

(One  Cluster  Approximation') 

Internal  Energy: 

The  internal  energy  of  the  Lenz-Ising  model  is,  in 
the  first  approximation, 

U  =  <H>  =  -MoBEi  <Pi  >  “  i  j  (  n.  n.  )  <PiPj  > 

=  -PoBNx  -  ^JNZx2 


where 


x  =  xi 1  and  Z  =  Zi  (  2 )  is  the  number  of  nearest 

neighbors . 

y  =  pair  correlation  coefficent  =  x2  in  this  1-cluster 
approximation  for  independent  clusters. 

Entropy: 

The  entropy  is 

S  -  S~(D  =  EiS~i<l>  =  NS1' .  (  1 ) 

S  =  -NkB  [Js(l+x)ln**(l+x)  +  1-x)  ln^(  1-x)  ] 

Free  Energy: 

The  Helmholtz  free  energy  is 
F  =  U  -  TS 
F  =  -  MoBNx  -  %JNZx2 

+  NkBT[^(l+x)ln^(l+x)  +  ^(l-x)ln^(l-x)] 
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It  is  convenient  to  work  with  a  dimensionless  version  of  the 


free  energy: 

3  =  (F/NJ)  =  -  Br  x  -  )£Zx2  +  Tt  [%( 1+x) ln%( 1+x) 

+  J5(l-x)lnJS(l-x)]  (6) 

where  Br  =  PoB/J  and  Tr  =  kBT/J  are  dimension¬ 
less  versions  of  the  magnetic  field  and  the  temperature. 

The  values  of  3  given  by  Eq.  (6), as  a  function  of  x, 
are  shown  in  Fig.  7  for  Tr  =  3,6,9,12,24,36.  The  critical 
temperature  for  this  first  approximation  is  Tr,ci  which  is 
the  same  as  the  pair  coordination  number  in  this 
approximation,  (see  Eq.  (8)).  The  values  selected  in  this 
figure  are  multiples  of  Tr  ,  c .  The  external  magnetic  field 
Br  , is  zero  in  Fig.  7. 

The  units  used  in  all  figures  for  3 , Br  , and  Tr  are 
called  "reduced"  units  and  are  dimensionless.  They  are 
given  in  equation  (6). 

Similar  results  are  given  in  Fig.  8  when  the 
external  magnetic  field,  Br  ,  is  equal  to  10.  This  is  a  very 
strong  field  and  is  equal  to  about  107  gauss  for 
ferromagnetic  systems.  The  minimum  free  energy  here  is 
close  to  saturation  (x=l),  as  compared  to  Fig.  7,  where  Br 
was  zero. 

Minimization*. 

The  equilibrium  state  of  this  system  is  found  by 
differentiating  F,  or  5 ,  with  respect  to  x: 
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FREE  ENERGY  PHI  vs.  MAGNETIZATION  X,  FOR  FIXED  VALUES  OF 
AND  MAGNETIC  FIELD  Br 
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[x  eqn] 

0  =  (d3/dx)  =  -  Br  -  Zx  +  Tr  [Jgln^d+x) 
-  Jgln^(l-x)] 


Hence : 

Br  +  Zx  =  l$Trln[(l+x)/(l-x)] 

Considered  as  an  equation  determining  x  =  x(Tr  , Br  ) ,  this 
is  transcendental  ,  that  is, the  methods  of  algebra  cannot 
solve  this  equation  for  x(Tr,Br).  However,  an  explicit 
result  is  possible  for  Tr  =  Tr (Br ,x): 

2  ( Br  +  Zx) 

Tr  — -  (7) 

ln[ ( 1+x) / ( 1-x ) ] 

Tables  and  graphs  for  Tr  =  Tr  (Br  , x)  can  be  inverted  to 
give  x  =  x  ( Tr  ,  Br  )  . 

The  equilibrium  free  energy  in  this  first 
approximation  is  calculated  by  solving  Eq.  (7)  for  Tr 
and  x  when  Br  is  kept  at  a  constant  value.  These  values  are 
used  in  Eq.  (6)  to  calculate  the  minimum  free  energy. 

The  results  are  given  in  Fig.  9  for  several  different 
values  of  Br  . 

There  is  a  critical  temperature  in  this  case,  which 
we  can  determine  by  the  conditions  that 
B=0 

and 

xsO  &  x<>0 


Equation  (7)  becomes 

2Zx  2Zx 

Tr,  cl  = -  =  —  =  Z 

ln[ ( 1+x) /( 1-x ) ]  2x 


(8) 
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There  are  two  branches  for  x  =  x(Tr  ,0)  when  Tr  <  Tr,ci .  3 

is  a  minimum  on  one  of  these  and  is  a  maximum  on  the  other. 
We  use  d2S/dx2  to  determine  which  is  which: 

d2S/dx2  =  -Z  +  Tr{JS[l/(l+x)]  +  Hi  l/(l-x)]} 

=  "Z  +  ^Tr[l/(l+x)  +  l/(l-x)] 

See  Fig.  7.  We  observe  that  for  Br  =0,  the  branch 
(x  =  0)  is  always  a  solution  of  (d<5/dx)  =  0.  However, 
this  solution  makes  §  ( x;  Tr  ,  Br  =0  )  a  minimum  only  for 
Tr  >  Tr.cl  =  Z.  For  Tr  <  Tr.c1  ,  this  branch  makes 
3 ( x; Tr  , Br  =0 )  a  maximum  ( i . e .  ,  d2 3 /dx2  <  0 ) .  Thermostatic 
equilibrium  occurs  only  when  F  is  a  minimum  vs .  changes  in 
x.  So,  we  select  x-branches  accordingly. 

Magnetization  x  is  compared  in  Fig.  10,  to  Tr  for 
various  values  of  the  magnetic  field  Br  .  Complete 
saturation  is  indicated  by  ±1,  with  all  the  spins  up  or 
down.  A  change  of  phase  occurs  when  Br  =  0 , and  this  is 
shown  to  occur  in  the  figure  for  Tr.c1  =  12,  and  also 
derived  in  Eq.  (8).  The  curve  for  Br  =  0  approaches  the  x  = 
0  line  from  either  side,  and  almost  crosses  the  line,  before 
turning  and  becoming  parallel  with  this  line  for  values  of 
Tr  >12. 

The  entropy  for  one  particle  clusters  is  presented  in 
Fig.  11,  as  a  function  of  the  reduced  temperature,  for 
various  values  of  the  magnetic  field.  The  entropy  in  bits 
per  particle  is  calculated  by  dividing  the  entropy 
( joules/oK)  by  ln(2) . 
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Tr  {REDUCED  TEMPERATURE) 

Fig.  10.  MAGNETIZATION  X  vs,  TEMPERATURE  Tr  FOR  VARIOUS  VALUES  OF 

MAGNETIC  FIELD  Bf 
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Fig.  12.  SPECIFIC  HEAT  vs.  REDUCED  TEMPERATURE  FOR  AN  FCC  LATTICE 
(1-CLUSTER  APPROXIMATION)  (See  Bragg  -  Williams  Fig.  5). 


The  specific  heat  Cb  is  compared  in  Fig.  12  with  Tr 
and  different  magnetic  fields.  These  results  can  be 
compared  directly  with  Fig.  5  from  reference  (17).  The 
results  in  the  first  approximation,  using  the  Morita 
expansion  of  the  free  energy,  and  truncating  the  series 
after  the  first  term,  are  the  same  as  the  Weiss  calculation 
for  one  particle  systems  using  self-consistent  field  theory. 
In  Fig.  5  (Ch.  II)  these  results  are  labeled  as  Bragg  - 
Williams,  which  is  a  refinement  of  Weiss's  calculations. 

A  contour  plot  of  the  entropy  in  "bits”  is  presented 
in  Fig.  13.  This  is  only  the  entropy  part  of  the  free 
energy  in  Eq.  (6)  and  all  values  are  positive  and  above 
the  x,y  plane  as  they  should  be.  This  plot  of  the  entropy 
will  be  referred  to  in  subsequent  approximations. 

Second  Order  Approximation:  Bethe-Peirels  Theory 

(Two  -  Cluster  Approximation) 

Internal  Energy: 

The  internal  energy  of  the  Lenz-Ising  model  is,  in 
the  second  approximation, 

U  =  <H>  =  -  PoBNx  -  ^JNZyi  (9) 

where  x  =  xi  (  1 )  ,  yp  =  xi  ,  i  +  p  <  2 )  , 

and  the  interaction  energy  is  represented  by  feiJNZyi . 

Compare  to  U  in  the  first  approximation  where  y  =  x2 . 

Entropy : 

The  entropy  is 

S  =  S*< 1 >  +  S‘<  2> 
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Fig.  13.  ENTROPY  FOR  1-PARTICLE  CLUSTERS  (UNITS  ARE 

BITS)  vs.  x,y. 
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i.e.  for  the  first  two  terms  of  the  free  energy  expansion 


where 

S‘d  >  =  Ei  =iMSi  *<  i  >  =  NS*d  > 

S*d)  =  -  NkB  [%( 1+x)  ln(  .  )  +  1-x)  ln(  .  )  ] 

S*<2)  =  E<  i  ,  j)  [Si  j<  2)  -  Sid)  -Sjd>] 

S*(2)  =  ^NkBSpZp  [S(  2)  i  ,  i  +P  -  2S(  1 )  i  ] 

S'd)  is  all  the  entropy  for  1-particle  clusters  and 
S~  (  2 )  is  all  the  entropy  for  2-particle  clusters, 
where  p  =  spacing  of  a  pair  of  particles 
Zp  =  pair  coordination  number 

=  number  of  pairs  of  spacing  p. 

S‘(2)  =  -^NkBE  p  Zp { [%( l+2x+yP ) ln( . ) 

+2(H) ( 1-yp ) ln( . )  +  *( l-2x+yP ) ln( . ) ] 
-2[*(l+x)ln(.  )  +  *(l-x)ln(  .  )]} 

Free  Energy: 

The  Helmholtz  free  energy  is 
F  =  U  -  TS 

5  =  F/NJ  =  (U  -  TS) /NJ 

3  =  -Bex  -  ^Zyi  +  Tr  {J$(l+x)ln(  .  )  +  ^(l-x)ln(.) 

+  5^2  P  Zp  {  [ ^ ( l+2x+yp  )  ln(  .  )  +  2(*) (1-yp )ln( .  ) 

+  l-2x+yP ) ln( . ) ] 

-  2[*S(l+x)ln(  .  )  +  $(l-x)ln(  .  )]}} 

The  equilibrium  state  of  this  system  is  found  by 
dif f erentating  §  with  respect  to  x  and  yi  : 

[x  eqn] 

0  =  (<S5/Sx) 

0  =  -Br  +  Tr  {%ln%(  1+x)  -  J$ln%(l-x) 


(10) 
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+  >2pZp{[2(34)ln%(l  +  2x+yp  ) 

-  2 (M) ln%( l-2x+yp ) ] 

-  2[^ln^(l+x)  -  Mln^(l-x)]}. 


1+x  l+2x+yp  1+x 

2Br/Tr  =  In - +  3£PZp[ln - -  2 In - ] 

1-x  l-2x+yp  1-x 


[yp  eqn] 

0  =  (>*5/<Fyp  ) 

0  =  -  %Z6'p\  +  ^Tr  {Zp  [^4ln^(  l+2x+yp  ) 

-^ln^(l-yp)  +  %ln^ln( l-2x+yp ) ] } 

Hence : 


( l+2x+yp ) ( l-2x+yp ) 

4Zi<5'pi  -  Tr  Zp  ln[ - ] 

(l-yp  )2 

We  must  solve  these  two  simultaneous  equations. 

Let  p>l.  Then  the  [yp  equation]  is  (for  Tr >0 ) 

[ ( l+2x+yP ) ( l-2x+yp ) 

0  =  In  - 

(1-7P  )2] 

since  Zp  =/0  for  any  p.  Hence, 

( l+2x+yp ) ( l-2x+yp ) 

-  =  1 

(l-yp )2 

(l+2x+yp ) (l-2x+yp )  =  (l-yp)2 

(l+2x)(l-2x)  +  [ ( l+2x)+( l-2x) ]yp  +  yp2  =  l-2yp+yp2 
4x2+2yp  =  -2yp 
from  which  we  find 

yp  =  x2  ,  p>l 

That  is,  in  the  "2-cluster  approximation",  spin-correlation 
extends  exactly  as  far  as  spin-interaction,  which  is 
"nearest  -  neighbor"  in  this  case. 

We  can  now  simplify  the  [x]  equation: 
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l+2x+yp  1+x 

E  P>  1  Zp  [In - -  21n - ] 

l-2x+yp  1-x 


l+2x+x2  1+x 

=  Ep>iZp[ln - -  21n ] 

1 - 2x+x2  1-x 

1+x  1+x 

=  E  P> 1 Zp [In  ( - ) 2  -  21n  — ]  -  0 

1-x  1-X 

Hence,  the  equations  become: 

1+x  l+2x+yi  1+x 

[x]  2Br/Tr  =  In - +  ^Zi  [In - 21n  ---] 

1-x  l-2x+yi  1-x 

( l+2x+yi ) ( l-2x+yi ) 

[y]  4/Tr  =  ln[ - ]  (11) 

( 1-yi  )2 

These  are  transendental  when  regarded  as  determining  x  &  y 
as  functions  of  Tr  &  Br  . 

Certain  results  can  be  extracted  from  [x]  &  [yi  ]  : 

1.  When  Br  =0,  inspection  shows  that  (x  =  0)  is  always 
a  solution  of  [x] .  We  then  use  [yi ]  to  determine 
yi  =  yi  ( Tr  ,  0  )  : 


(1+71  )2 


4/Tr  =  In  - 

(1-yi  )2 

1+yi 

implies  2/Tr  =  In  -  =  2arctanh  yi 

l-yi 

yi  =  tanh  1/Tr 

See  the  following  derivation  for  this  result. 
1+y  1+y 

[Let  a  =  In  -  ,  -  =  exp(a)  =  b  ,  1+y  =  b-by  , 

1-y  1-y 


b-1  exp(a)-l 


by+y  =  b-1  ,  (l+b)y  =  b-1  ,  y  =  -  =  - 

b+1  exp(a)+l 
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exp( a/2 ) -exp( a/2 ) 

-  =  tanh(a/2) 

exp ( a/2 ) +exp ( a/2 ) 

a  =  2arctanh  y 

See  figure  14  and  15  for  a  plot,  of  x  vs .  Tr  and  y  vs. 

Tr  . 


An  alternative  method  was  actually  used  to  generate 
the  results  plotted  in  figure  14  and  15.  This  is  the 
Simplex  method  of  minimization,  and  does  not  use  the  methods 
of  calculus  to  achieve  a  minimum  value  for  a  function[6] .  In 
this  approximation  it  is  used  as  a  convenient  method  to 
obtain  the  results  shown  in  these  figures.  The  same  results 
would  be  obtained  by  solving  the  minimum  equations  for  [x] 
and  [y]  given  by  Eq.  (11).  It  will  be  shown  in  later 
approximations  that  it  is  necessary  to  use  the  Simplex 
method  to  mimimize  the  free  energy  equations  that  are 
obtained  by  truncating  the  Morita  expansion  of  the  free 
energy.  It  will  be  shown  in  these  higher  approximations 
that  the  mimima  occur  on  the  boundary  of  the  function,  and 
minimization  using  the  methods  of  calculus  will  not  work  in 
these  situations.  The  Simplex  program  for  this  second 
approximation  is  listed  in  Appendix  E. 

The  equilibrium  free  energy  is  given  in  Fig.  16  as  a 
function  of  Tr  for  various  values  of  Br  .  If  this  result  is 
compared  with  Fig.  9  in  the  1-cluster  result  it  will  show 
that  the  equilibrium  free  energy  results  are  the  same  in 
both  cases,  at  least  to  the  accuracy  of  the  grid  size  used 
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Fig.  15.  PAIR  CORRELATION  MOMENT  y  vs.  TEMPERATURE  Tr  FOR  VARIOUS  VALUES 
OF  MAGNETIC  FIELD  Br  (1,2-CLUSTER  APPROXIMATION) 


to  plot  these  results.  See  Appendix  G,  Table  2  for  the 
actual  values  of  the  equilibrium  free  energy. 

2.  There  is  a  critical  temperature.  We  can  deterrmine  it 
as  we  did  for  the  1-cluster  approximation,  by  setting 
the  conditions 

Br  =  0  and  (  x~  0&x=/0  ) 

[x]  4x 

0  =  2x  +  ^Zi  [ - -  4x] 

1+yi 

1 

implies  yi  =  - 

(Zi-1) 

as  can  be  shown  below. 


1 

0  =  1  +  Zi  [ - -  1] 

1+yi 

1  1  ( 1+yi  )  -  1  yi 

Zi  1+yi  1+yi  1+yi 

yi  +1  1 

Zi  = - =  1  +  — 

yi  yi 


l  1 

Zi  -  1  =  —  ,  yi  = - 

yi  Zi  -1 

That  is,  at  Tr,cii  ,  yi  =  ( Zi  —  1 ) - 1  .  We  use  this  in 

[yi ]  to  determine  Tr  ,  c  : 

2  2  2 

Tr  ,  cl  I  = - = - - - - - 

1+yi  1  1  Zl  ~1  +  1 

In -  ln[  ( 1+  — )(1 - )-!]  In - 

1-yi  Zl  -1  Zl  -1  Zl-1-1 


2  2 


Zl  1 

In -  In - 

Zl  -  2  1  -  2  ( Zi  )-l 
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Tr  (REDUCED  TEMPERATURE) 

Fig.  16.  EQUILIBRIUM  FREE  ENERGY  $  vs.  TEMPERATURE  Tr  FOR  VARIOUS 
VALUES  OF  MAGNETIC  FIELD  Br  (1,2-CLUSTER  APPROXIMATION) 


Expanding  this  as  a  series  in  1/Zl  allows  convenient 
comparison  with  Tr.c1  '• 


Tr , ci  I 


2 

-2/Zi  -  y2{  2/Zi)2  -  1/3  ( 2/Zi  )3  +... 


2 

(  -2/Zi  )  [  1  +Jg(2/Zi)  +1/3  (  2/Zi  ) 2  +...] 

Zi 

1  +%( 2/Zi  )  +1/3  ( 2/Zi  )2  +.  .  . 

=  Zi{l  -  [%( 2/Zi  )  +1/3  (  2/Zi  ) 2  +...] 

+  [%( 2/Zi  )  +. .  -]2  +  -  -  ■} 

=  Zl{l  -$(2/Zi)  +(-1/3  +1/ 4 )  ( 2/Zi  )2  +...} 

=  Zl{l  -1/Zl  -(1/3)  (1/Zl  )2 
=  Zi  -l-(l/3)  (1/Zl  )  +.  .  . 

=  Tr  ,  cl  -1  -d/3)  (1/Tr  ,  CI  )  +.  .  . 

=  Tr,CI[l  -1/Zl  -(1/3)  (l/(Zl  )2  -...] 

As  in  the  "1-cluster"  approximation,  the  critical 
temperature  depends  on  the  crystal  lattice  structure 
including  its  dimension  --  only  through  Zi  =  the  number  of 
nearest  neighbors,  12  (See  Fig.  1,  Ch.  I).  So,  for  example, 
the  two  dimensional  hexagonal  lattice  has  the  same  behavior 
as  the  three  dimensional  cubic,  in  this  approximation. 

Since  exact  calculations  show  that  their  behaviors  are 
different,  we  note  that  a  dependence  on  only  Zi  is  a  failing 
of  this  approximation. 

It  is  well-know  that  these  sorts  of  approximations 
become  more  exact  as  Zi  — >  *  (i.e. ,  as  dimension  >  ® ) , 
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and  becomes  very  bad  as  Zl  — >2  (i.e.,  as  dimension  >1)- 

We  now  examine  the  stability  of  the  solutions  : 

>5'2  5  Tr  1  1  1  1 

3xx  = - =  --  [ - + - +  EpZp  (  + 

£x2  2  1+x  1-x  l  +  2x+yp  l-2x+yP 

1  1 

- )] 

1+X  1-X 


S'  2§  Tr  11 

3  7  P  x  = - -  — SpZp  (  -  ) 

6'yptSx  4  l+2x+yp  l-2x+yp 

5  25  Tr  1  1 

§xyp  - - =  -~Zp  (  -  ) 

SxSy  p  4  l+2x+yp  l-2x+yp 

S'  25  Tr  1  2  1 

5  7  p  y  p  = - -  — Zp  (  +  +  ) 


S'  yp  2  8  l  +  2x+yp  1-yp  l-2x+yp 

Evaluating  these  on  the  branch  '• 


gives : 


Br  =  0 

x  =  0  &  yi  =  tanh  1/Tr 

yp>i  =0 


$25  Tr  2 

---  =  —  [2  +  Zi  ( - -  2)  +2  p  >  1  Zp  ( 0  )  ] 

S  x2  2  1+yi 


1 

Tr  [1  +  Zl  ( - 1)] 

1+yi 


S'  25  Tr 

- =  —  Zp  ( 0 )  =  0 

SxSyp  4 

$ 25  Tr  Zp  1  ( 1+yi  )_  1  +  1  ( 1-yl )_  1  p=1 

£yp2  4  0  P>1 


Inspection  shows  that  the  branch  Br =0  yp =tanh  1/Tr  p-1 

x=0  0  p>l 

is  stable  (  3 ’ ’ >0  )  for  Tr  >  Tr.cll  ,  but  is  unstable 
(§xx  <  0  )  for  Tr  <  Tr  ,  c*  I  . 

There  is  another  branch  for  Br  -  0  &  x  = /  0  ; 
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it  is  stable  for  Tr  <  Tr , c1 1  .  We  now  study  it  .  .  . 
Equation  [x]  can  be  solved  for  yi  =  yi(x;Br/Tr)  : 


l+2x+yi 

1+x 

2Br 

1+x 

[In - 

-  21n 

— ] 

-  ( - - 

In  — )(Zx/2) 

l-2x+yi 

1-X 

Tr 

1-x 

l+2x+yi 

1 

2Br 

1+X 

1+x 

In 

- - 

2[  —  ( 

- - 

In  — ) 

+  In - ] 

l-2x+yi 

Zi 

Tr 

1-x 

1-x 

l+2x+yi 

2Br 

1 

1+x 

In 

- - 

2  [ - 

+  (1 

- )  In 

—  3  =  a 

l-2x+yx 

Zl  Tr 

Zi 

1-x 

- >  ( l+2x+yi )/( l-2x+yi  )  =  expa  =  b 

- >  l+2x+yi  =  b-2bx+byi 


yi -byi  =  b-2bx-l-2x 
(l-b)yx  =  (b-l)-2x( 1+b) 

- >  yi  =  -l-2x[ ( l+b)/( 1-b) ]  =  -1  +2x[( 1+expa) /( -1+expa ) ] 

yi  =  -1  +2x[ (expa+1 )/(expa-l ) ] 
yi  =  -1  +  2x/(tanh{2Br /Zi  Tr  +  ( 1-1/Zi  )  ln[  ( l+x)/(  1-x)  ]  }  )-i 

Suppose  Br  =  0  and  Zi - >  ®  Then 

yi  =  -1  +  2x/(tanh{ln[ ( 1+x) /( 1-x) ] } )_  1 

=  -1  +2x/[4x/( 2+2x2 ) ]  =-l  +  1  +  x2  =  x2 
as  we  would  expect. 

Suppose  Br  =  0  .  Then 

yi  =  -1  +  2x/( tanh{ ( 1-1/Zl ) ln[ ( 1+x) /( 1-x ) ] } )_ 1 
Also  for  x  ~  0,  we  find: 

yi  =  -1  +  2x/tanh[  ( 1-1/Zl  )  (  2x)  ] 

=  -1  +  2x/( 1-1/Zl )2x 

=  -1  +  1/(1-Zi)  =  (-1  +  1/Zi  +  1)/(1  -  1/Zi  ) 

=  l/(Zi  -  1)  ,  as  expected. 

These  results  are  plotted  in  Fig.  17,  with  yi  vs  x. 
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We  observe  a  deviation  from  the  1-cluster  result,  y  -  x2 ,  as 
a  result  of  including  pair-correlation  effects. 

This  result  for  yi  =  yi  (x;Br/Tr  )  can  be  used  in  the 
[yi ] -eqn , 

[yi  ]  ( l+2x+yi  )  ( l-2x+yi  ) 

4/Tr  =  In - 

(1-71  )2 

to  calculate  Tr  as  a  function  of  x  &  Br/Tr.  The  results, 

yi  =  yi  (x;  Br  /Tr  ) 

Tr  =  Tr  (x;  Br/Tr  ) 
can  be  inverted  to  give 

x  =  x  ( Tr  ,  Br  ) 
yi  =  yi  (Tr  ,Br  )  . 

of  course,  7p>i  =  x2  . 

[Results  ...  in  the  same  pattern  as  for  1-cluster  case.] 
Although  the  results  in  this  1-  and  2-cluster 
approximation  seem  to  be  reasonable  e.g.  a  real  critical 
temperature  is  calculated,  it  was  discovered  that  already  at 
this  level  of  approximation  there  is  a  problem  with 
truncating  the  free  energy  expansion.  The  results  of  this 
difficulty  will  show  up  very  clearly  in  the  next  higher 
cluster,  i.e.  the  3-cluster  approximation,  where  for 
example,  a  complex  value  is  obtained  for  the  critical 
temperature.  The  following  figures  will  illustrate  the 
roots  of  the  problem. 

In  Fig.  18  we  plot  only  the  internal  energy  part  of 
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Fig.  17.  1,2-CLUSTER  APPROXIMATION  FOR  y  COMPARED  TO  1  CLUSTER 


Eq.  10  for  the  free  energy,  which  is  the  same  as  given  by 

Eq.  9.  The  contour  plotting  method  used  in  these  results  is 

similar  to  those  discussed  in  [23] ,  [24] , [25]  .  In  Fig.  19 
only  the  entropy  part  of  the  free  energy  is  plotted.  This 
shows  some  very  unphysical  behavior  where  the  entropy 
becomes  negative  for  some  regions  of  the  x,y  plane.  This  of 
course  is  physically  impossible,  and  will  be  a  particular 
problem  in  the  next  higher  approximations . 

Figures  20  and  21  show  the  entire  free  energy  given 

by  Eq.  10.  It  shows  two  minima,  and  the  illustrates  the 

distortion  in  the  free  energy  surface  due  to  the  behavior  of 
the  entropy.  This  distortion  will  become  worse  with  the 
third  approximation,  and  the  minima  will  occur  on  the 
boundary,  where  the  methods  of  calculus  cannot  be  used  to 
find  a  minimum.  This  2-cluster  approximation  has 
illustrated  the  effect  of  truncating  the  free  energy 
expansion,  although  the  difficulties  resulting  from  this 
truncation  are  more  evident  in  the  3-cluster  results.  A 
listing  of  the  program  HIDDEN6  used  in  this  contour  plotting 
is  given  in  Fig.  33  of  Appendix  E. 

Third  Order  Approximation:  Limited  3-Cluster 

We  will  restrict  ourselves  to  a  "square"  lattice  and 
include  only  nearest  neighbor  and  next-nearest  neighbor 
pairs  and  retain  only  the  most  compact  triplet  and  ignore 
all  other  n-clusters  (n  1  2). 
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Fig.  19.  ENTROPY  (IN  BITS)  FOR  1,2-PARTICLE  CLUSTERS  vs.  x,y 


Fig.  20.  5  (INTERNAL  ENERGY  &■  ENTROPY)  FOR  1,2-PARTICLE 
CLUSTERS  vs.  x,y  FOR  Br=0,  Tr=10.  (  $  CLIPPED  AT  1  >  2) 
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In  this  restricted  view  we  have  the  following 
"  square"  lattices  *. 

In  one  dimension,  the  lattice  is  linear,  and  the 
retained  pairs  are  given  by  p=l  and  p=2.  The  quantities 
p,q,r  are  the  distances  between  particles.  The  retained 
triangle  is  described  by  p-1,  q^l)  r“2.  The  coordination 
numbers  in  1-dimension  are: 

pairs - >  Zp  ;  Z(  2)  1  =2,  Z(  2 )  2  =2 

compact  triplet - >  Zpqr  ;  Z( 3 ) x 1 2  =  6. 

In  two  dimensions,  the  lattice  is  square  and  the 
retained  pairs  are  p=l,and  p=2  (or  strictely,  p=^ 2  but  we 
have  retained  p  as  an  "index"  of  distance).  The  retained 
triangle  is  given  by  p=l,  q=l ,  r=2  ( actually  r=-4*2.  The 
coordination  numbers  are: 

pairs - >  Zp  ;  Z(  2)  1  =  4,  Z(  2>  2  =  4 

triplet - >  Zpqr;  Z(3)ii2  =  24 

In  three  dimensions  the  lattice  is  cubic.  The 
retained  clusters  are  like  those  in  the  square  case,  and  the 
coordination  numbers  are: 

jpairs - >  Zp  ;  Z(2)i  -  6,  Z(  2 )  2  -  12 

triplet - >  Zpqr;  Z(  2  )  1 1  2  -  72 

These  results  can  also  be  obtained  arithmetically,  using  the 
rule  that  particles  are  located  at  lattice  sites  according 
to 

r->ii,  i2.  .  .  id  =  aZ  Dotsi  niocb~o< 
where  the  b"«x  are  orthogonal  unit  vectors  for  "square 


98 


lattices 


a  is  the  lattice  constant,  and  the  nicx  are 


integers . 


We  select  the 

partial 

it  ii 

e 

at  r-  > 

A 

o 

1! 

and  vary 

it  to  locate  other  particles 

The  distance  between  " . " 

and 

any  other  particle  is 

;  j ln| <  2)  =2 

i  =  1  D  ni  2 

We  calculate 

the 

coordination  numbers 

of  each 

particular 

sort  of 

n-cluster  by 

counting  the  number  of  In’s 

that 

generate  n-clusters  of 

that 

particular  sort. 

In  one  dimension  ,  D= 

1 ,  and  we  have : 

{n=+l ,  n=-l}  <--> 

smallest  2-cluster 

<  —  > 

Z(  2)  1 

=  2 

{n=+2 ,  n=-2}  <  — > 

next 

smallest  2-cluster 

< 

:  —  >  Z(  2)  x  = 

2 

n  n*  p 

q 

r 

1  : 

2  1  . 

2 

1 

2 

L  2 

1 

1 

-1  +1  1 

1 

2 

Z<  3  )  1 1  2 

=  6 

+1 

1  1 

1 

2 

-l 

2  1 

2 

1 

-2 

1  2 

1 

1 

In  two  dimensions,  D 
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and  we 
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nx 

ny  1 

Ml 

0 

1 

i 

0 

-1 

i 
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=  4 

1 

0 

i 

-1 

0 

i 

1 

1 

42 
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We  continue  in  a  similar  way  as  in  the  first  two 
approximations . 

Internal  Energy: 

The  internal  energy  is, 

U  =  <H>  =  -MoBNx  -  ^JNZiyi 

Entropy: 

The  entropy  is, 

S  =  S*( 1 >  +  S‘( 2)  •  +  S‘( 3>  1 

(where  the  prime  notation  reminds  us  that  we  are  restricting 
our  attention  to  n.n.  pairs  and  the  n.n.n.  pairs  and  only  to 
triangles  for  the  triplets.) 

The  intrinsic  entropy  for  1-clusters  is, 

S^i)  =  EiS*(i)i  =  NS*  (  1 )  =  -Nkb  {*£(  l±x)  In.  } 

(The  notation  "In."  means  that,  the  quantity  that 
immediately  precedes  the  function  "In"  is  repeated  as  its 
argument . ) 

For  pairs  the  intrinsic  entropy  is, 

S~  (  2 )  ’  =  S*  (  2 )  .  ,  .  +1  +  S*  <*).,.  +2 

=  -NkB{Zi {[(1/4) (l+2x+yi ) In . +etc] -2 [ ] } 

+  Z2 { [ ( 1/4 ) ( l+2x+y2 ) In. +etc] -2 [ ] }  } 

For  triangles  the  intrinsic  entropy  is, 

S~<  3) ’  =  S*( 3)  . , .+1, .+2 

=  -NkB  Zi  i  2  {  [  ( 1/8)  ( l+3x+2yi +y2+zi  l  2  )  In.  +etc]-etc  } 
(Note:  In  this  description  the  used  as  a  subscript, 

as  in  S'*.  ,  is  meant  to  replace  the  i  usually  used  to 

represent  a  single  particle.  Also  in  the  pair  notation  the 
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i,j  would  be  replaced  by  .,.+1  and  n.n.n.  pairs  are  written 
as  . , .+2.  For  triangles  we  use  . , . +1, .+2  instead  of  the 
usual  i , j , k. ) 

Free  Energy  •• 

The  Helmholtz  free  energy  for  this  limited  3-cluster 
approximation  is 

S  =  (F/NJ)  =  [-Br  x-l/2Zl  Yl  ]  +  Tr  {[J*(l+x)ln.  +  (k(l-x)ln.] 

+%£p  =  i ,  2  Zp  {  [%( l  +  2x+yp  )  In.  +  2*^4  ( 1-yp  )  In . 

+^(l-2x+yP  ) In .  ]  -2 [J£(  1+x) In.  +HK  1-x) In .  ]} 

+  ( 1/6 )Zn  2  {[( 1/8)  ( l+3x+2yi  +y2  +zl  l  2  )  In. 

+2*( 1/8 ) ( l+x-y2 -zi  l  2 ) In . 

+  ( 1/8 ) ( l-x-2yi +y2  +zi  l  2 ) In . 

+  ( 1/8  )  ( l+x-2yi  +y2  -  zi  l  2  )  In . 

+2* (1/8) (l-x-y2+zii2 ) In . 

+  ( 1/8 ) ( l-3x+2yi  +y2 -zi  l  2 ) In. ] 

-2[54(  l+2x+yi  )  In.  +2*?4(  1-yi  )  In.  l-2x+yi  )ln.  ] 

-  [34(l+2x+y2  )  In. +2*^(  l-y2  )  In .  +%( l-2x+y2  )ln.  ] 

+3C^(l+x)ln.+^(l-x)]>} 

Minimization : 

We  will  minimize  the  free  energy  using  calculus 
minima  and  also  with  the  Simplex  algorithm,  reference  (6). 

The  equilibrium  state  of  this  system  is  determined  by 
finding  the  minimum  values  for  x,yi,y2,zii2  which  are 
obtained  by  calculating  the  derivative  of  $  first  wrp.  to  x, 
and  then  yi , etc . ,  setting  the  eqns .  equal  to  zero,  and 
solving  for  x ,  yi  ,  y2  ,  zi  l  2  . 

Calculating  the  derivative  of  §  wrp.  x,  and  setting 
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"the  result  equal  to  zero,  is  called  the  [x]  eqn. 


[x]  o  =  («/tfx) 

l+x 

—  —  Br  +  Tr  {*£ln 

1-x 


l+2x+yp  l+x 

+S(Ip  =  i,  2  )ZP[*ln - -  In  ---  ] 

l-2x+yp  1-x 

l  +  3x+2yi +y2 +zi  i  2  l+x-y2-zii2 

+( 1/6 ) Zi l 2  [ ( 3/8 ) In  -  +  (2/8) In  - 

l-3x+2yi +y2 -zi  l  2  l-x-y2+zii2 

( l+x-2yi +y2 -zi  l  2  )  l+2x+yi 

( 1/8 )  In - -  In - 

( l-x-2yi +y2 +zi  l  2  )  l-2x+yi 

l+2x+y2  l+x 

-  ^ln - +  (  3/2 )  In - ]} 

l-2x+y2  1-x 

l+x  l+2x+yp  l+x 

2*Br/Tr  =  In - +  Jg(Zp  =  i,  2  )Zp[ln - -  2 In - ] 

1-x  l-2x+yp  1-x 

Zi  l  2  l  +  3x+2yi  +y2  +zi  l  2  l+x-y2-zii2 

+ - [31n -  +  21n - 

24  .  l-3x+2yi +y2 -zi  1 2  l-x-y2+zii2 

l+x-2yi +y2 -zi  l  2  l+2x+yi 

+  In -  -  81n - 

l-x-2yi +y2 +zn  2  l-2x+yi 

l+2x+y2  l+x 

-  41n -  +  121n - ] 

l-2x+y2  1-x 


[yi] 

0  =  {Si/Syi  ) 

=  -4Zi  +  Tr  (Zi  [ln^(l+2x+yi  )  -  21n%(l-yi)  +  ln%( l-2x+yi  ) ] 
+  2(4/6)Zii2{[(2/8)ln(l/8)( l  +  3x+2yi  +y2  +zi  l  2 ) 

-  2(1/8) In (1/8) ( l-x-2yi  +y2  +zi  l  2 ) 

-  2(1/8) In (1/8) ( l+x-2yi  +y2-zii2  ) 

+  (2/8)ln(l/8) (l-3x+2yi  +y2 -zi  l  2 )] 
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-  %[ln%( l  +  2x+yi  )  -  21n^(l-yi)  +  ln%( l-2x+yi  ) ] } } 
( l+2x+yl  ) ( l-2x+yi  ) 


4/Tr  =  In  - 

(l-yi  )2 

Zi  l  2  ( l+3x+2yi  +y2  +  zi  l  2  )  ( l-3x+2yi  +y2  -zi  l  2  ) 

+ - [In - 

3Zi  ( l-x-2yi  +y2  +zi  l  2  )  ( l+x-2yi  +y2  -zi  l  2  ) 

( l+2x+yi ) ( l-2x+yi ) 

-  21n - ] 

( l-y2 )2 


[y2  ] 

0  =  (S$/Sy  2) 

=  Tr  {^Z2  [^ln  %( l+2x+y2 )-2*^ln%( l-y2 )+%ln%( l-2x+y2 ) ] 
+  ( 1/6  )  Zi  l  2  [  ( 1/8  )  ln(  1/8  )  ( l+3x+2yi  +y2  +zi  l  2  ) 

-2(1/8) In (1/8) ( l+x-y2 -zi l 2 ) 

+  ( 1/8) ln( 1/8 ) ( l-x-2yi +y2  +  zi  l  2 ) 

+( 1/8 ) ln( 1/8 ) ( l+x-2yi +y2 -zi l 2 ) 

-2 ( 1/8 ) ln( 1/8 ) ( l-x-y2  +zi l 2 ) 

+(1/8 ) ln( 1/8 ) ( l-3x+2yi +y2 -zi 12)] 

- [^ln^( l+2x+y2 )-2^1n^(l-y2 )+^ln^( l-2x+y2 )]} 

( l+2x+y2 ) ( l-2x+y2 ) 

—  >  0=(J$)($)Z2ln - 

(1-72 ) 2 

Zi  l  2  ( l+3x+2yi  +y2  +zi  x  2  )  ( l-x-2yi  +y2  +zi  1  2  ) 

+ - [In - 

6*8  ( l+x-y2 -zi  1  2  )2  ( l-x-y2 +zi  x  2  )2 

( l+x-2yi  +y2  -zi  1  2  )  ( l-3x+2yi  +y2  -zi  1  2  ) 

*  - - - 

( l+x-y2  -zi  1  2  )2  ( l-x-y2  +zi  12)2 

( l+2x+y2 ) ( l-2x+y2 ) 

-21n  -  ] 

(1-72 ) 2 

( l+2x+y2 ) ( l-2x+y2 ) 

-->  0  =  In  - 

(1-72 ) 2 
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Zi  1  2  ( l  +  3x+2yi  +y2  +zi  l  2  )  ( l-3x+2yi  +y2  -zi  l  2  ) 

+ - [In - 

6Z2  (  l+x-y2 -ZI  1  2  )2  (  l-X-y2 +Z1  1  2  )2 

( l-x-2yi  +y2  +zi 1  2  )  ( l+x-2yi  +y2  -zi  l  2  ) 

* - 

(  l+X-y2  -zi  12)2  (  l-x-y2  +  Z1  12)2 

( l+2x+y2 ) ( l-2x+y2 ) 

-21n  -  ] 

( l-y2 )2 


[zi  1  2  ] 

0=  ($a/j  zn2) 

=  Tt  {  ( 1/6  )  Zi  i  2  [  ( 1/8  )  ln(  1/8  )  ( l+3x+2yi  +y2  +zi  l  2  ) 
-2 ( 1/8 ) ln( 1/8 ) ( l+x-y2 -zi  l  2 ) 

+  (1/8)  In (1/8)  ( l-x-2yi  +y2+zii2  ) 

-  ( 1/8 )  ln(  1/8  )  ( l+x-2yi  +y2-zii2  ) 

+2(1/8) ln( 1/8 ) ( l-x-y2  +zi  l  2 ) 

-  ( 1/8 )  ln(  1/8 )  ( l-3x+2yi  +y2-zii2  ]} 

( l+3x+2yi  +y2  +zi  l  2  )  ( l-x-2yi  +y2  +zi  l  2  )  ( l-x-y2  +zi  l  2  )2 

1=  - 

( l-3x+2yi  +y2  -  zi  l  2  )  (l+x-2yi  +y2-zii2  )  (l+x-y2-zii2  )2 


We  begin  the  study  of  these  by  supposing  that  Br  =  0. 
By  inspection,  we  find  that  (  x=0,zii2=0  )  then  satisfy  [x] , 
[ zi  l  2 ]  identically: 

1+yp 

[x]  0  =  lnl  +  ^(Zp  =  i,  2 )Zp [In - -21nl] 

1+yp 


Zi  1 2  l+2yi  +y2  l-y2  l-2yi+y2 

+ - [31n - +21n - In - 

24  l+2yi  +y2  l-y2  l-2yi+y2 

1+yi  l+y2 

-81n - -41n - +61nl] 

1+yi  l+y2 


[zi  l  2  ] 
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( l  +  2yl  +72 ) ( l-2yi  +72 ) ( I-72 )2 

- -  1 

( 1+271 +72 ) ( 1-271 +72 ) ( 1-72 )2 

Equations  [71]  &  [72]  become: 


[71]  4  I+71  Zi  1 2  (I+271+72)  I+71 

-  =  2  In - +  - [  21n -  -41n - ] 

Tr  I-71  3Zl  ( 1-271 +72  )  l-yi 

2  I+71  Zl  1 2  I+271+72  I+71 

-  =  In - + - [In - -21n - ] 

Tr  1-71  3Zl  1-271 +72  I-71 

[y2  ] 

1+72  Zl  1  2  (  1  +  271  +72  )(  1-271  +72  )  1+72 

0  =  21n - + - [21n  - -41n - ] 

I-72  6Z2  (1-72)2  I-72 


1+72  Zi  1 2  ( 1+271 +72 )( 1-271 +72 )  I+72 

0  =  ln - + - [In - -21n - ] 

I-72  6Z2  (1-72)2  I-72 

We  will  now  solve  [72]  for  71  =  71  (72  ;  Zi  1  2 /Z2  )  ,  and 


use  the  result  in  [71]  to  obtain  Tr  =  Tr  (72  ;  Zi  1  2 /Zi  ;  and 
Zi  1  2  /Z2  )  . 

[y2  ] 

6Z2  1+72  ( 1+271 +72 )( 1-271 +72 )  1+72 

---  [-In - ]  =  In - ' -  -  21n - 

Z112  I-72  (1-72)2  I-72 

( l+2yi +72  )  ( I-271 +72  )  6Z2  I+72 

ln - =  [2 - ]  [In - ] 

( I-72  ) 2  Zi  1 2  I-72 

( l+2yi  +72  )( I-271  +y2  )  6Z2  I+72 

- =  exp  { [2 - ]  [  In - ]}  =  b 

( 1-72)2  Z112  I-72 

( l+2yi +72  )  ( I-271 +72  )  =  (l~72)2b 

(I+271  )  (I-271  )  +  [(l+2yi  )  +  (l-2yi  )]y2  +  y22  =  bd-72)2 
1  -  471 2  +  272  +  72 2  =b( I-272  +  72 2  ) 

-4yi2  =  [l-2y2+y22]b  -  1  -  2y2  -  y22 
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yi  2  =  ^[1  +  2y2  +  y2  2  -  b  +2by2  -by2  2  ] 
yi  2  =  %[(l-b)  +  2  ( 1+b  )y2  +  (l-b)y22] 
yi  2  =  14 [  1  — b 3  [1+72  2  +2  [  ( 1+b)  /( 1-b)  ]y2  ] 

We  have  now  derived  the  equations  for  b  and  yi 2 . 

To  evaluate  b  and  yi 2  we  have  the  following  calculations 
in  one  and  two  dimensions. 

In  One  Dimension  Z2  =  2  and  Zi  1 2  =  6 .  Hence 

b  =  exp{ [2-6*2/6] ln[ ( l+y2 )/( l-y2 ) ] } 

=  exp{[  0  ] 

=  1  and 

yx  2  =  ( 1/4 ) [  0  +  2*2y2  +  0*72  2 ] 
yi  2  =  y2 

Summarizing:  b  =  1 

y2  =  yi 2 

In  Two  Dimensions  Z2  =  12  and  Zi  1 2  =  24.  Hence 
b  =  exp{ [2-6*12/24] ln[(l+y2 )/(l-y2 )  } 

=  exp{ [2-3] ln[( I+72 )/( I-72 )] 

=  exp{-ln[ ( I+72 ) /( I-72 ) ] } 

=  [(1+72  )/(l“72  )]  (-D 
=  [(1-72 )/( 1+72 )] 

and  the  relation  between  71  and  72  is, 
yi  2  =  %[(l-b)  +  (l-b)y22  +2(l+b)y2  ] 

=  ^{(l-[(l-72 )/(l+72 )] )  +2( l+[ ( I-72 )/ ( I+72 ) ] )72 
+( l-[ ( 1-72 )/( 1+72 ) ] )  } 

=  (272  )/( I+72)  +  2  (  2/(  I+72  )  ) 72 

+( 272 /( 1+72 ) )722  } 
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=  %[2y2 /( l+y2 ) ] [3+72 2 ] 

=  *z[y2 /( l+y2  )  ]  [3+y22  ] 

Summarizing:  b  =  [ ( l-y2 ) / ( l+y2 ) ] 

yi2  =  ^[72/(1+72  )]  [3+y22  ] 
We  use  these  results  in  [yi  ]  to  obtain  Tr  . 


[yi] 

2 

1+yi  Zi  l  2 

l+2yi  +y2 

1+yi 

— 

=  In - + - 

[In - 

-  21n  ---] 

Tr 

1-yi  3Zi 

l-2yi  +y2 

1-yi 

We  derived  the  result  for  yi 2  as: 

yi 2  =  ( l-b)+2( l+b)y2 +( l-b)y22 } 

6Z2  l+y2 

where  b  =  exp{[2-  - ] In  - }  (from  p.105) 

Zi  1  2  l-y2 

Using  these  equations  we  found  in  ONE  dimension  that 
y2 =yi 2  . 

In  ONE  dimension,  if  we  use  this  result  in  [yi  ]  above 
and  recall  that  Zl  1 2  =  6 , Zi  =  2  then  we  derive  yi  as: 

2  1+yi  6  l+2yi +yi 2  1+yi 

-  =  In - +  -  [In - -  21n - ] 

Tr  1-yi  6  l-2yi+yi2  1-yi 

1+yi  1+yi  1+yi 

=  In - +  21n - -  21n - 

1-yt  1-yi  1-yi 

1+yi 

=  In - +  0 

1-yi 

-->  2/Tr  =  ln[  ( 1+yi  )  /( 1-yi  )  =  a 

-  >  ( 1+yi  )  /  ( 1-yi  )  =  exp  (a)  =  b 

-->  1+yi  =  b  -  byi 

-  -  >  ( 1  +b )  yi  =  b  -  1 
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—  > 


71 


—  > 

—  > 

—  > 


71  2 


b-1  exp(a)-l 
b+1  exp(a)+l 
exp(a/2) -exp( -a/2 ) 
exp ( a/2 ) +exp ( -a/2 ) 
=  tanh(a/2) 


=  tanh(l/Tr  ) 

We  know  that  this  is  the  same  result  as  in  the 
2-cluster  case. 

( 1/4 ) [ - (  ) 2y2  +  2(2)72  ~(  )722] 

=  (1/4) [4(l  +  l/2(  ) 72  -(  ) 72 2 ] 

since  (  )y 2 2  =0 


712  =  ( 1/4 ) [4 ( 1+1/2 (  ) 72 ] 

=  [1+1/2 ( 2-6 Z2 /Zl 12 ) ]72 
=  [2  -3Z2/ZH2  ]72  =  [2-(3*2)/6]72 
In  one  dimension,  we  found  712  —  72 ,  and  we  find  it  again 
from  this  result. 

Hence,  in  one  dimension,  71  =  tanh(l/Tr  )  . . .same  as  in  the 

2-cluster  approximation. 


Summar7 

For  an7  dimension,  we  find,  for  Br =0 

6Z2  I+72 

b  =  exp{  [2 - ]  [In - ]} 

Z112  I-72 


1+b 

71  =  J^{(l-b)  [l+2(---)72+722  ]} 

1-b 


2  I+71  Zl  1 2  I+271 +72  I+71 

-  =  In - + - [In - -  21n - ] 

Tr  I-71  3Zi  1-271+72  1-71 


x  =  zi  1 2  =  0 
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The  existence  of  negative  values  for  the  entropy  in 
the  truncated  free  energy  was  illustrated  in  the  previous 
section,  see  Fig.  19.  The  effect  on  the  free  energy  surface 
was  shown  in  Figs.  20  and  21.  In  the  third  approximation 
these  negative  values  of  the  entropy  have  severe  effects  on 
the  derived  results.  If  a  critical  temperature  is  derived 
from  these  minimized  equations  it  will  have  a  complex  value 
[23].  To  examine  the  entropy  values  resulting  from 
truncating  the  expansion  after  the  third  term  the  values  are 
shown  in  Figs.  22  and  23  for  two  different  values  of  z,  the 
triple  correlation  moment  .  In  Fig.  22  the  value  of  the 
entropy  values  are  given  as  a  function  of  magnetization,  x, 
and  the  pair  correlation,  y,  for  a  constant  value  of  z  = 

0.2.  Similar  results  are  given  in  Fig.  23  for  z  =  0.  The 
entropy  values  shown  in  the  upper  half  of  both  figures  are 
negative  in  some  cases  and  greater  than  unity  for  some 
values.  Neither  of  these  kind  numerical  are  physically 
realistic:  the  entropy  must  be  il  and  can  never  be  negative 
for  this  system.  The  entropy  values  where  obtained  with  the 
program  GRID  6  and  the  listing  is  given  in  Appendix  E. 

An  attempt  at  reducing  the  magnitude  of  these  non¬ 
physical  results  is  shown  in  the  lower  half  of  figures  22 
and  23.  This  reduction  is  suggested  by  numerical  techniques 
which  only  use  %  of  the  value  of  the  last  term  in  a  series 
in  order  to  speed-up  convergence.  We  have  used  only  ^  of 
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x  (MAGNETIZATION) 


a  (TRIPLE  CORRELATION  MOMENT) 


1 

.8 


-.6 

-.8 


-1 


3-2 


-.42 

.6 

.61 

.09 

-.87 

.4 

.79 

.64 

.40 

.09 

-.57 

.80 

.81 

.71 

.61 

.52 

.45  “ 

.72 

.83 

.80 

.75 

.72 

.75 

.851.06 

.2 

.63.79 

.80 

.76 

.73 

.73 

.80 

.971.271.79- 

.53 

.75 .78 

.74 

.68 

.64 

.63 

.69 

.821.06 

0 

.39 

.70 

.77.73 

.64 

.55 

.46 

.41 

.39 

.38  - 

-.32  .56 

.76 

.76.66 

.51 

.35 

.19 

.03 

-.17 

-.2 

.77 

.74.58 

.35 

.09- 

-.21 

-.58 

— 

.50 

.73.51 

.16 

-.2G  • 

.78 

-.4 

— 

.57.49- 

-03 

-.72 

— 

.49— .17 


y  (PAIR  CORRELATION  MOMENT) 


Z1 11=48  - 


1 


.8 


3=2 


.6 


.4 


.2 

0 

-.2 

-.4 

-.6 


-.8 


-1 


0 

.2 

.4 

.6  .8  1 

.09 

- 

.66 

.31 

-.33 

.79 

.65 

.38 

.02 -.60 

.83 

.79 

.63 

.40 

.13-21 

.83 

.85 

.75 

.59 

.39 

.15-11  -.42 

.80 

.87 

.83 

.71 

.54 

.34 

.12  -.13  -  40  -.68— 

.73 

.87 

.87 

.79 

65 

.47 

.26 

.02  -24  -.55 

.83 

.88 

.84 

.73 

.58 

.38 

.14 

-  14-.48 

.85 

.87 

.81 

.67 

.49 

.25 

-04 

-.41 

.82 

.85 

.77 

.60 

.37 

.08 

-.31 

— 

.62 

.81 

.72 

.52 

.23 

-17 

.64 

.67 

.42 

.02 

— 

.55 

.31 

INCLUDES  ONLY 
Vi  THE  3-CLUSTER 
ENTROPY 


Fig.  22.  ENTROPY  FOR  1,2,3-PARTICLE  CLUSTERS  vs.  x,y  WITH^=.2 
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x  (MAGNETIZATION 


3  (TRIPLE  CORRELATION  MOMENT) 


1  r- 
.8  - 
.6 

.4 
.2 
0 

-.2 
-.4 
-.6 

-.8 

-1 


3=0 


.39 

-.97 

.76  .37 

-.25 

-1.33 

.85  .77  .49 

.14 

-  .30 

-1.13 

.82 

.92  .81  .62 

.43 

.22 

-  .02 

.72 

.96 

.95  .85  .73 

.63 

.54 

.50 

.47 

.52 

.91 

.99 

.96  .89  .81 

.75 

.73 

.77 

.90  1.12 

-.10  .69 

.94  1.00 

.97  .90  .83 

.78 

.79 

.86 

1.03  1.34  1.87 
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the  entropy  contributon  from  the  three  particle  cluster  (the 
last  custer  term  in  the  truncated  series)  and  re-calculated 
the  total  entropy  for  the  1,2,3-particle  clusters.  The 
results  do  show  an  improvement  in  the  entropy  values',  some 
become  less  negative  and  all  are  i 1 .  In  Fig.  23  the  entropy 
value  is  unity  for  z-y=x=0,  as  it  should  be  for  both 
calculations . 

The  equilibrium  free  energy  is  plotted  in  Fig.  24 
for  the  third  approximation  as  a  function  of  Tr  .  The 
values  were  determined  by  the  Simplex  minimization 
algorithm.  A  listing  of  the  Simplex  program  is  give  in 
Appendix  E,  and  the  essential  steps  are  given  in  reference 
[6], 

As  mentioned  in  the  previous  section  (4.1)  the 
minimum  occur  on  the  boundary  and  the  calculus  minimazation 
will  not  work.  We  have  used  the  Simplex  algorithm  to 
minimize  the  truncated  free  energy  equation.  This 
minimization  technique  is  not  a  calculus  type  minimization 
method  and  can  be  used  to  obtain  minima  which  are  located  on 
the  boundary. 

In  Fig.  24  the  equilibrium  values  are  denoted  by  "X” 
for  5  when  no  correction  is  made  to  the  entropy  contribution 
of  the  third  term  in  the  truncated  free  energy.  These 
values  for  §  are  very  different  than  the  solid  curves  which 
are  the  equilibrium  values  for  the  free  energy  that  were 
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Fig.  24.  EQUILIBRIUM  FREE  ENERGY  *  vs.  TEMPERATURE  Tr  FOR  VARIOUS 
VALUES  OF  MAGNETIC  FIELD  Br  (1,2,3-CLUSTER  APPROXIMATION) 


obtained  in  the  earlier  approximations  and  re-plotted  in 
Fig.  24.  It  was  noted  previously  that  the  numerical  value 
of  3eq.  would  not  change  very  much  for  the  higher  order 
approximations  as  compared  to  the  approximations  for  1 12- 
particle  clusters.  To  evaluate  the  effect  of  only  using  % 
of  the  entropy  contribution  of  the  third  term,  as  was  done 
in  Figs.  22  and  23,  the  values  of  3  were  re-calculated  using 
only  ^  of  the  third  term.  These  results  are  denoted  by  the 
large  black  dots,  and  show  good  agreement  with  the  earlier 
results . 


Fourth  Order  Approximation: ( 4-Cluster ) 

The  following  is  a  list  of  the  principal  steps  in  this 
derivation . 

The  Helmholtz  Free  Energy  is  given  as  F  =  U  -  TS. 

1.  Calculate  U. 

2.  Calculate  S.  ~  indicates  intrinsic  entropy. 

s  =  S*(1)  +  S“<2)  +  S*<3)  +  S“<«>  (11) 

S‘<1)  is  the  same  as  in  the  1-c  approximation. 

S“(2)  "  2-c 

S~<3)  "  3-c 

S"(4)  =  S(4)  -  4S(3)  +  6S(2)  -  4S(1) 

S(  1 )  ,S(2)  ,S(3)  are  known  from  the  earlier 

approximations  and  are  given  as  R1,R2,R3,  respectively  in 

the  program  SIM4,see  the  listing  in  Appendix  E. 

S(  4)  is  determined  from  P(  4 )  (JJ  ’  .  .  .M  ’  ’  ’  )  which  is  the 

probability  distribution  for  only  4  particles.  (See 

section  3.3).  The  equations  for  this  distribution  are 
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denoted  by  R4  in  SIM4 ,  (see  Appendix  E) .  (Note  that  "R"  is 
used  to  represent  entropy  in  SIM4. )  The  total  entropy,  R,  is 
given  as  follows 

R  =  R1  +  [Zl/2] (R2-2R1 )  +  [Zlll/6] (R3-3R2+3R1 ) 

+  [Zllllll/24] (R4-4R3+6R2-4R1)  (12) 

Z1  is  the  pair  coordination  number 
Zlll  is  the  triplet  coordination  number 
Zllllll  is  the  4-particle  coordination  number 

Z1,Z111  are  known  from  the  previous  approximations.  Zllllll 

is  the  only  new  coordination  number. 

3.  The  Helmholtz  free  energy  is  given  as, 

5  =  -Brx  -  %Zi  y  +  Tr  (R)  (13) 

5,  Br  and  Tr  are  the  reduced  free  energy,  magnetic  field 
and  temperature  respectively. 

4.  Minimize  the  free  energy,  and  obtain  the  minimum  value  of 

5 .  The  probabililty  distribution  associated  with  this 
minimum  value  of  3  is  the  equilibrium  macrostate  for  this 
system.  From  this  macrostate  all  quantities  of 
thermodynamic  interest  can  be  derived  by  the  standard 
formulas . 

The  following  is  a  description  of  the  4-cluster 
approximation  results  for  the  FCC  lattice. 

Internal  Energy: 

U  =  <  Hamiltonian  >  =  -  poBNx  -  MJNZ(  1  )  <P  l  P  2  > 

These  quantities  have  been  previously  defined  and  the 
definitions  are  repeated  here. 

po  is  the  magnetic  moment  of  each  particle 


115 


B  is  the  external  magnetic  field 
N  is  the  number  of  particles  in  the  system 
x  is  the  magnetization  per  particle 
The  interaction  energy  is  given  by  the  second  term. 

J  is  the  exchange  integral 

Z1  is  the  pair  coordination  number 

Pi  is  the  microstate  of  one  of  the  particles  in  the 
pair  and  P2  is  the  microstate  of  the  other 
particle.  This  pair  of  particles  are  nearest 
neighbors  (n.n. )  and  the  average  value  of  the 
product  of  the  two  microstates  is  called  yi  ,  the 
pair  correlation  coefficient. 

U  can  be  written  as , 

U  =  N [  -P o  Bx  -  J^JZiyi] 

Entropy  '• 

The  entropy  for  this  4-cluster  approximation  is  given 

by, 

s  =  S~<1>  +  S‘(2)  +  S‘(3)  +  S‘(4) 

The  intrinsic  entropy  for  each  cluster  is  represented  by  S 
with  a  superscript  "  ~  ",  and  each  cluster  approximation  by 
a  number,  i.e.  S' ( 2 )  is  the  intrinsic  entropy  for  the  second 
approximation.  The  intrinsic  entropy  for  each  approximation 
is  written  as  follows. 

S~ < 1 )  =  S(  1 )  for  the  first  approximation,  since  there 

is  only  one  particle  in  each  cluster  and  it  is  not  necessary 
to  correct  for  the  prescence  of  any  other  clusters.  S(  1 )  is 
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given  by  the  entropy  for  the  one-particle  exact  calculation, 
see  section  3.0. 

S*(i)  =  SO)  =  -kB[J$(l+x)ln*(l+x)  +  *S(l-x)lnJg(l-x)]  (14) 

The  intrinsic  entropy  for  the  2-cluster  is  written  as, 

S'  <  2  )  =  £  i  ,  j  [S<  2  )  i  ,  j  -  SO  )  i  -  SO  )  j  ] 

S' (2)  =  J^NIpZp  [S(  2)  i ,  i  +p  -  2SO)i] 

S(  2 ) i ( i +  p  is  the  entropy  derived  in  the  exact  calculation 
for  two  particles,  see  section  3.1,  and  is  given  as, 

SO)  =  -kB  [%( l+2x+yp  )  ln(  .  )  +  J£(  1-yp  )  ln(  .  )  +  l-2x+yp  )  ln(  .  )  ] 

(15) 

SO)  is  the  same  as  above.  Then  S'O)  is  given  as, 

S'  (  2  )  =  -^NkB2PZp{[)4(l+2x+yp  )ln(  .  )  +  %(l-yp)ln(.) 

+  %( l-2x+yp ) ln( . )  ]  -  2[^( 1+x) ln( . )  +  %( 1-x) ln( . ) ] }  . 

The  intrinsic  entropy  for  the  3-cluster  is  written  as, 
S'  (  2 )  ’  =  S<3)  -  2S(2)p  =  i  -  S(  2)  p  =  2  +  3SO) 

S< 3 >  is  the  entropy  derived  in  the  exact  calculation  for 
three  particles,  see  section  3.2,  and  is  given  as, 

S(3)  =  -kB  [  ( 1/8 )  ( l+3x+2yi +y2 +zi  l  2  )  ln(  .  ) 

+  2( 1/8 ) ( l+x-y2 -zi l 2 ) ln( . ) 

+  ( 1/8) ( l-x-2yi +y2 +zi  l  2 ) ln( . )  +  ( 1/8 ) ( l+x-2yi +y2 -zi l 2 ) In ( . ) 

+  2 ( 1/8 ) ( l-x-y2 +zi  i  2 ) ln( . )  +  ( 1/8 ) ( l-3x+2yi +y2 -zi  l  2 ) ln( . ) ] 

(16) 

S(  2 )  and  S(  1 )  are  the  same  as  given  above.  The  intrinsic 
S'O)  is  given  in  the  previous  section. 

The  intrinsic  entropy  for  the  4-cluster  is 
S'O)  =  SO)  -  4S<3)  +  6S(2)  -  4SO> 

All  terms  in  this  equation  have  been  previously  defined. 
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SM)  ,  the  exact  distribution  for  only  4  particles,  is  given 
in  section  3.3  S(  1  )  ,  S(  2 )  f  S(  3  )  are  given  by  Eqs  . 

( 14 ) , ( 15 ) , ( 16 )  respectively. 

The  total  entropy  for  the  first,  second,  third  and 
fourth  cluster  is  given  by  Eqs.  (11), (12).  In  SIM4  the 
total  entropy  for  these  four  clusters  is  denoted  by  R. 

Free  Energy: 

The  free  energy  is  given  by  Eq.  (13)  and  is 
denoted  by  line  190  in  SIM4  (Appendix  E,  Fig.  32). 

Minimization : 

The  minimization  of  the  free  energy  is  calculated 
with  SIM4  and  the  results  presented  in  Fig.  26. 

Discussion: 

The  entropy  values  for  the  first  four  clusters  are 
shown  in  Fig.  25  in  a  similar  way  as  the  entropy  values  are 
given  in  Figs.  22  and  23  for  the  third  order  approximation. 
(These  values  were  calculated  with  the  program  GRID6  that  is 
listed  in  Appendix  E  as  Fig.  34).  The  entropy  as  a  function 
of  x  and  y  is  given  in  Fig.  25  for  z  =  w  =  0  in  the  upper 
half  of  the  figure,  and  for  z  =  w  =  0 . 2  in  the  lower  half. 

In  this  case  however  it  is  not  necessary  to  adjust  these 
entropy  values  as  was  done  in  the  third  order  approximation. 
The  values  are  plausible  and  physically  acceptable  -  none 
are  >1  and  none  are  negative.  Note  also  that  for  x  =  y  =  z 
=  w  =  0  the  entropy  has  a  maximum  value  of  unity  as  it  should. 
The  equilibrium  free  energy  is  given  in  Fig.  26  as  a 
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Fig.  26.  EQUILIBRIUM  FREE  ENERGY  $  vs.  TEMPERATURE  Tr  (4-CLUSTER) 


function  of  the  reduced  temperature  Tr  and  is  denoted  by 
"x".  The  reduced  magnetic  field  is  zero.  The  solid  curve 
is  the  equilibrium  free  denergy  for  a  1-particle  clusters 
and  as  previously  noted  these  equilibrium  vaues  will  not 
change  very  much  with  higher  order  approximations.  The 
equilibrium  values  are  also  given  in  Table  1,  Appendix  G. 

For  any  value  of  Tr  ,  the  minimization  program  SIM4  may 
calculate  values  for  equilibrium  §  which  are  more  negative 
than  those  plotted  in  Fig.  26  and  represented  by  "x".  The 
correlation  values  x,y,z,w  for  these  more  negative  i3  values 
will  be  less  acceptable  than  the  correlation  values 
associated  with  the  less  negative  value  of  § .  Table  2  in 
Appendix  G  gives  values  of  x,y,z,w  for  Tr  =  4,8,10,16.  The 
values  of  §  plotted  in  Fig.  26  also  agree  very  closely  with 
the  equilibrium  5  values  for  the  lower  order  approximations. 
The  reason  that  other  values  are  calculated  for  3  is  that 
these  other  equilibrium  values  are  associated  with  the 
spurious  minimum  that  are  located  on  the  boundary.  These 
values  of  5  are  too  low  and  the  values  of  x,y,z,w  associated 
with  them  are  physically  unacceptable.  The  entropy 
calculated  for  these  $  values  will  be  negative.  These  were 
the  criteria  used  in  selecting  the  §  values  that  are  plotted 
in  Fig.  26. 


121 


APPENDIX  A 


COORDINATION  NUMBERS 


Linear  Lattice:  Isingl.BAS 
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Square  Lattice: 
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Cubic  Lattice: 
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"4-cube"  Lattice: 
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5-cube"  Lattice: 
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"6 -cube"  Lattice: 
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APPENDIX  B 


ISING1  PROGRAM 

Regarding  the  program  ISING1.BAS  for  computing 
Z(2)p--pair  coordination  numbers,  we  list  the  following 
rules.  (See  Fig.  27) 

1.  Put  one  particle  at  the  origin  of  the  lattice.  Put 
the  other  particle  at  the  location  specified  by  the  lattice 
numbers  N1,N2, . . . , ND  =  |N  (where  D  =  dimension  of  the 
lattice) . 

Since  the  lattice  is  generated  by  the  basis  vectors 
Bl_> ,  B2->  , . . . ,BD_>  =  ! B“ >  by  the  rule 

r~>  -  Z  i  Nl  Bl  ■  >  , 

then  the  location  (of  the  particle)  specified  by  j  N  has 
cartesian  coordinates 

r  =  E  i  Ni  ( Bi  )  -  >  j 

where  (Bi  )->j  are  the  cartesian  coordinates  of  |B-> 

This  is  the  location  rule. 

The  "n-cube"  lattices  are: 


n 

1 

2 

3 

4 

5 

6 


n  -  cube 
linear  lattice 
square  lattice 
cubic  lattice 
hyper  -  cube  lattice 
who  knows  what  it  is  called? 
etc . 
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These  lattices  are  generated  by  the  "standard” 
orthonormal  basis,  whose  cartesian  components  are 
(Bl  )j->  =  tfji  . 

The  distance  rule  is  given  by  the  following.  The 
distance  between  the  particle  at  the  origin,  and  the  one 
specified  by  |N  is 

P  =  !  !r->-  0->  !  |  =  !  |r->  !  1 

=  £j  [EiNi  (Bi  ->  )j  ]2 

For  "n-cube"  lattices,  this  is 

P  =  -J  [S  j  (Si  Nitfji  )2] 

=  -J-  [2  j  Nj  2  ]  . 

2.  Step  over  all  values  of  |N, 

reject  P  if  (a)  P  =  0  (we  have  accidentally  put  the 

roving  particle  at  the.  origin) 

(b)  P  >  Pmax  (we  have  an  jN  that  puts  the 
roving  particle  too  far  from  the 
origin)  . 

3.  We  keep  a  list,  P(I),  of  distinct  values  of  P,  always 

arranged  in  ascending  order.  We  keep  a  second  list,  C(I), 
of  the  number  of  times  the  distance  P(I)  has  been  obtained. 
These  lists  are  initialized  to  {  C(I)  =  0  } 

{  P(I)  =  Pmax+1  }  . 

Upon  moving  the  roving  particle  to  each  new  location, 
and  computing  P,  we  scan  the  list  P(I)  starting  from  the 
smallest  value  (at  the  top).  For  each  value  of  I  (1=0  to 
Imax),  we  compare  P  with  P(I): 
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(a)  if  P  =  P(I),  then  we  increment  the  counter  C(I), 
C(I)+1,  and  then  quit  the  comparison. 

(b)  if  P  <  P(I),  then  we  recognize  that  we  have 
encountered  a  new  value  of  P,  which  we  must  insert  into  the 
list.  So  we  pause  to  open  up  a  space  in  the  list,  by  moving 
all  values  from  the  current  value  P(I)  to  the  top  of  the 
list  P(Imax)  down  by  one.  And  we  also  move  the  counters. 

(Of  course,  P(Imax)  and  (Imax)  fall  off  the  end  of  the 
lists . ) 


For  J  =  Imax  to  I  step-1 
P(J+1)  =  P(J) 

C ( J+l )  =  C(J)  . 

Next  J 

We  then  insert  P  at  location  I ,  and  set  its  counter  to  1 : 
P(I)  =  P 
C(I)  =  1  . 

We  then  quit  the  comparison. 

(c)  if  P  does  not  match  any  value  on  the  list,  we  ignore 

it  . 

4.  After  the  roving  particle  has  moved  over  its  entire 

domain,  we  print  the  results: 

For  I  =  0  to  Imax 

Print  I,  P(I),C(I) 

Next  I 

Of  course,  we  can  read  this  as 
p,  Z(  2  )  p  . 
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150  REM  N=2 
160  REM 

170  OPTION  BASE  1 
172  CLEAR 

174  POKE  18,0: POKE  19,0: POKE  20,0: 

180  DIM  B ( 3 , 3 ) 

200  REM  DEFINITION  OF  BASIC  VECTORS 
210  REM  B(I,X)  is  Xth  component  of  Bi 
220  B( 1 , 1 ) =1  :  B( 1 , 2 ) =0  :  B(l,3)=0 
230  B( 2 , 1 ) =0  :  B(2,2)=l  :  B(2,3)=0 
240  B( 3 , 1 ) =0  :  B(3,2)=0  :  B(3,3)=l 
250  MAX=2:M2=MAX*MAX 
260  DIM  R( 3 ) , C(M2 ) 

300  FOR  I=-MAX  TO  MAX 
310  FOR  J=-MAX  TO  MAX 
320  FOR  K=-MAX  TO  MAX 
330  FOR  L=1  TO  3 

340  R(L)=I*B(1,L)+J*B(2,L)+K*B(3,L) 
350  NEXT  L 

360  D=R ( 1 ) *R ( 1 ) +R ( 2 ) *R ( 2 ) +R ( 3 ) *R ( 3 ) 

370  IF  (D=0)  OR  (D>M2 )  THEN  400 

380  C(D)=C(D)+1 

400  NEXT  K 

410  NEXT  J 

412  PRINT  ”I  =  "  ;  I 

420  NEXT  I 

490  OPEN  #1,"P:M  OUTPUT 
500  PRINT  #1, "D  " ; MZ(D) M 
510  FOR  D=1  TO  M2 
520  PRINT  #1,D;"  M;C(D) 

522  PRINT  D; ”  " ;C(D) 

530  NEXT  D 
600  CLOSE  #1 


FIGURE  27.  Program  Isingl  for  Calculating  Coordination 
Numbers 
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APPENDIX  C 


ISING2  PROGRAM 

Regarding  the  program  ISING2.BAS  for  computing 

Z(3)pqr  - triangle  coordination  number.  Put  one  particle 

oc  at  the  origin  of  the  lattice.  Put  another  ft  at 
a  location  specified  by  |M  and  a  third  r  at  a 
location  specified  by  |N.  Then 

roc-  >  =  0-> 

rtf  -  >  =  Z  i  Ni  Bl  -  > 

rT  -  >  =  Z  i  Nl  Bl  "  > 

The  distances  between  the  particles  are 
Docd  =  !  !  r<x_  >  -  rtf->  j  \  =  \  ;rtf->  \  j 
=  siti  cum  (Bi->  )j  ]2> 

Dar  =  !|ra->  -  rT~>  j  |  =  ||rT->!J 

=  J{Zj  [ZiNl  (Bl  - >  )i  ]2} 

Qir  =  ! !rtf->  -  rT->  j  \ 

=  J{Zj  [Zi  (Ml  -  Ni  )  (Bl -->  )j  ]2  } 

For  "n-cube"  lattices  these  are: 

Dafl  =  -J  (S  jMj2  ) 

Dar  =  -J(ZjNj2  ) 

Dtfa  -  J[Z  j  (Mj  -  Nj  ) 2  ] 

2.  Step  0>  and  r  throughout  the  entire  region 
surronding  a.  Arrange  Daft ,  Dar , and  Dftr 
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into  ascending  order: 


(Docfl  ,Dxr  ,Dflr  )  -->  (P,Q,R) 

i.e.  ascending  order  where  P  <  =  Q  <  =  R. 
At  each  pair  of  locations,  compute  Dxf'J  ,  Dar,Dftr. 


Go  immediately  to  the  next  location  if 

(a)  any  of  these  are  0  (at  least  one  particle  has 
accidently  been  placed  on  top  of  another) . 

(b)  any  of  these  exceed  Dmax  (at  least  one  particle 
has  gotten  too  far  away). 

3.  There  are  only  certain  possible  distances  between 
lattice  points.  We  know  what  these  are  from  the  "Z2p" 
work.  Call  these  possible  distances,  A,B,C, . . . ,  and 
arrange  them  in  ascending  order.  Then  the  combinatorially 
possible  3-tuples  are 


AAA 

AAB 

ABB 

AAC 

ABC 

ACC 

AAD 

ABD 

ACD 

AAE 

ABE 

ACE 

- 

- 

- 

- 

- 

- 

- 

- 

— 

BBB 

BBC 

BCC 

BBD 

BCD 

BDD 

BBE 

BCE 

BDE 

- 

- 

- 

- 

- 

- 

- 

- 

- 

ccc 

CCD 

CDD 

CCE 

CDE 

CEE 

- 

- 

- 

_ 

_ 

— 

Of  course,  not  all  of  these  combinatorially  possible 
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p-tuples  is  a  possible  triangle  on  a  given  lattice;  for 
example,  in  1-dimension,  Dsfi  +  DGr  =  Dsr .  In  practice  we 
are  only  interested  in  "small"  triangles.  So  in  practice, 
we  will  cut  off  each  sequence  when  the  corresponding 
triangle  is  "too  big" .  We  can  now  assign  each  remaining 
3-tuple  a  counting  number  I,  and  a  counter  C(I).  We 
increment  this  counter  whenever  P,Q,R  passes  the  appropriate 
entrance  requirement. 

For  n-cube  lattices,  the  possible  distances  are  m, 
m=l,2,3,...  .  (Not  all  values  of  m  occur  for  smaller  n-cube 

lattices . ) 

So  ■ 

FOR  J  =  0  TO  4 

D(J)  =  SQR(J) 

NEXT  J 

loads  a  more-than-adequate  set  into  D( . ) .  We  might  prefer 
to  do  this  "by  hand"  for  each  explicit  lattice,  to  ensure 
that  only  distances  possible  for  that  lattice  occur. 

Anyway,  D( . )  is  to  contain  the  possible  distances. 

Here  is  a  possible  scheme,  for  n-cubic  lattices,  that  cuts 
off. 
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TABLE  1 


TRIANGLE  COORDINATION  NUMBER 


1  1 

1  p  1 

[  9  1 

R 

0 

DCO) 

DCO) 

0(0} 

1 

0(0) 

0(0) 

0(1) 

2 

0(0) 

0(0) 

0(2) 

3 

D(Q) 

0(0) 

0(3) 

4 

DCO)  i 

I  0(1)  I 

f  0(1) 

5 

0(0)  1 

0(1) 

i  D(  2 ) 

6 

0(0)  | 

l  0(0) 

1  D{3) 

7  1 

1  0(0)  1 

r  D{2>  | 

r  0(2) 

8 

|  0(0)  . 

■  0(3) 

9 

|  D(Q)  | 

j  0(3) 

|  0(3) 

10 

[  D(l)  1 

1  0(1) 

r  o(D 

11 

,  Dd)  ; 

1  0(1) 

j  0(2) 

12 

1  ! 

1  0(2) 

[  0(2) 

13 

[  0(2) 

|  0(2) 

\  D(2) 

See  figure  23* 
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100  REM  THIS  PROGRAM  COMPUTES 

110  REM  N-COORDINATION  NUMBERS  (FCC) 

130  REM  4  DEC  1985  (copied  from  ISING3.msb  16  feb  1985) 

140  REM 
150  REM  N=3 
160  REM 

170  OPTION  BASE  1 
172  CLEAR 
174  TIME=0 
176  F$="scrn: " 

178  OPEN  F$  FOR  OUTPUT  AS  #1 
195  DEFINT  A-Z 

200  REM  DEFINITION  OF  BASIS  VECTORS 

210  REM  B(I,X)  is  the  Xth  component  of  Bi 

212  DIM  B( 3 , 3 ) 

214  REM  FACE  CENTERED  CUBIC 

220  B( 1 , 1 ) =1  :  B( 1 , 2 ) =1  :  B(l,3)=0 

230  B( 2 , 1 ) =0  :  B(2,2)=l  :  B(2,3)=l 

240  B( 3 , 1 ) =1  :  B(3,2)=0  :  B(3,3)=l 

250  MAX=2 

260  M2=MAX*MAX 

270  DIM  C (M2 , M2 ) 

280  DIM  R1  (3) ,R2(3) 

300  FOR  I1=-MAX  TO  MAX 
302  FOR  J1=-MAX  TO  MAX 
304  FOR  K1=-MAX  TO  MAX 
306  FOR  L=1  TO  3 

308  R1  (L)=I1*B(1,L)+J1*B(2,L)+K1*B(3,L) 

310  NEXT  L 

312  D1=R1 ( 1 ) *R1 ( 1 ) +R1 ( 2 ) *R1 ( 2 ) +R1 ( 3 ) *R1 ( 3 ) 

314  IF  (D1=0)  OR  (D1>M2 )  THEN  450 

320  FOR  I 2= -MAX  TO  MAX 

322  FOR  J2=-MAX  TO  MAX 

324  FOR  K2  =  -MAX  TO  MAX 

326  FOR  L=1  TO  3 

328  R2(L)=I2*B(1,L)+J2*B(2,L)+K2*B(3,L) 

330  NEXT  L 

332  D2=R2 ( 1 ) *R2 ( 1 ) +R2 ( 2 ) *R2 ( 2 ) +R2 ( 3 ) *R2 ( 3 ) 

334  IF  (D2=0)  OR  (D2>M2)  THEN  400 
338  D3=0 

340  FOR  L=1  TO  3 

342  R3=R1(L)-R2(L) :D3=D3+R3*R3 

344  NEXT  L 

FIGURE  28.  Program  Ising2  for  Calculating  Coordination 
Numbers 
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346  IF  (D3=0)  OR  (D3>M2)  THEN  400 

348  GOSOB  6000 

400  NEXT  K2 

402  NEXT  J2 

404  NEXT  12 

450  NEXT  K1 

460  NEXT  J1 

470  NEXT  II 

500  PRINT  #1, "Dl,  D2 ,  Z(D1,D2)" 

510  FOR  Dl=l  TO  M2 
512  FOR  D2=l  TO  M2 

522  PRINT  #1,D1; "  ";D2;M  ";C(D1,D2) 

524  PRINT  Dl;"  ";D2;"  ";C(D1,D2) 

530  NEXT  D2.D1 

534  PRINT  #1,  "This  calculation  took  "; TIMES 
536  PRINT  #1: PRINT  TIMES 
600  CLOSE  #1 
5000  P=D1 

5010  IF  D1>D2  THEN  T=D1 : D1=D2 : D2=T 
5020  IF  D2>D3  THEN  T=D2 : D2=D3 : D3=T 
5030  IF  D1>D2  THEN  T=D1 • D1=D2 : D2=T 
5040  C(D1,D2,D3)=C(D1,P2,D3)+1 
5050  D1=P 
5060  RETURN 

6000  D(1)=D1:D(2)=D2:D(3)=D3:P1=D1 
6010  M=3 

6020  FOR  G=M-1  TO  1  STEP  -1 
6030  FOR  F=1  TO  G 

6040  IF  D(F) >D(F+1 )  THEN  SWAP  D(F),D(F+1) 

6050  NEXT  F 
6060  NEXT  G 

6070  D1=D( 1 ) : D2=D( 2 ) : D3=D( 3 ) 

6080  C(D1,D2)=C(D1,D2)+1 
6085  D1=P1 
6090  RETURN 


FIGURE. 28  continued 
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APPENDIX  D 


CLUSTERS  OF  4-PARTICLES 

This  appendix  describes  "the  calculation  of  the 
coordination  number  for  a  cluster  of  four  particles.  Two 
pieces  of  information  are  obtained  from  the  calculation: 

(1)  the  shape  of  the  most  compact  figure  for  4-particles  and 

(2)  the  number  of  these  most  compact  figures  which  is  called 
the  cordination  number.  The  particles  are  fixed  in  position 
at  the  lattice  sites  of  the  crystal  (see  Fig.  1,  Ch.  I) 

Fig.  29  is  a  diagram  of  the  method  used  in  ISING4F 
to  calculate  the  distances  between  the  particles.  Use  one 
of  the  particles  as  the  origin  and  draw  a  vector  to  each  of 
the  other  three  particles  which  are  called  the  "rovers"  and 
labeled  1,2,3.  Calculate  the  distances  corresponding  to 
these  3  vectors  and  the  distances  between  each  of  the  3 
"rovers".  The  vectors  and  the  associated  distances  are 
labeled  as  shown  in  Fig.  29. 

The  calculation  of  the  six  distances  are  calculated 
and  stored  by  ISING4F  for  each  particle  that  is  choosen  as 
an  origin.  The  distances  are  sorted  and  the  number  of 
similar  figures  are  counted.  The  results  for  the  4-cluster 
case  are  given  in  Fig.  30  and  show  that  the  most  compact 
figure  is  an  equilateral  pyramid  with  six  sides  equal  and 
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that  there  are  48  of  them.  The  fact  that  all  six  sides  are 
equal  length  is  an  important  distinction.  This  has  been 
true  for  the  calculations  in  this  study.  It  will  not  be 
true  for  higher  order  clusters  e.g.  5-cluster  and  6-cluster 
configurations.  This  means  that  the  entropy  calculation 
will  be  done  in  a  different  way  than  described  in  this  work. 

Fig.  31  is  a  listing  of  the  steps  in  the  program 
ISING4F  to  carry  out  this  calculation  for  an  FCC  latice. 
Other  lattices  are  possible  by  changing  the  values  in  lines 
160  to  180  to  represent  other  lattice  sites. 
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4-CLUSTER 
3 "ROVERS” 


ORIGIN 


R1 

D1 

R2 

D2 

R3 

D3 

R?2 

D12 

RT3 

D13 

R23 

D23 

Fig.  29.  DIAGRAM  OF  A  4-PARTICLE  CLUSTER  SHOWING  THE 
ORIGIN  AND  3  "ROVERS”.  THE  VECTORS  R1  THRU  RZ3  AND 
CORRESPONDING  DISTANCES  D1  THRU  D23  ARE  USED  IN  THE 
PROGRAM  ISING  4F  TO  CALCULATE  THE  NUMBER  OF  THE  MOST 
COMPACT  FIGURES  FOR  AN  FCC  LATTICE. 
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PYRAMIDS  WERE  COUNTED. 


10  REM  THIS  PROGRAM  COMPUTES 
20  REM  N-COORDINATION  NUMBERS  (FCC) 

40  REM  4  DEC  1985  (converted  from  ISING3.msb  16  feb  1985) 
50  REM  see  pg.  18/VIII  and 
60  REM  N=4 
70  REM 

80  OPTION  BASE  1 
90  CLEAR 
100  TIME=0 
102  F$= " lptl : " 

104  OPEN  F$  FOR  OUTPUT  AS  #1 
110  DEFINT  A-Z 

120  REM  DEFINITION  OF  BASIS  VECTORS 

130  REM  B(I,X)  is  the  Xth  component  of  Bi 

140  DIM  B( 3 , 3 ) 

150  REM  FACE  CENTERED  CUBIC 

160  B( 1 , 1 ) =1  :  B( 1 , 2 ) =1  :  B(l,3)=0 

170  B( 2 , 1 ) =0  :  B( 2 , 2 ) =1  :  B(2,3)=l 

180  B(3,l)=l  :  B( 3 , 2 ) =0  :  B(3,3)=l 

190  MAX=2 

200  M2=MAX*MAX 

210  DIM  C (M2, M2, M2, M2, M2, M2) 

220  DIM  R1 ( 3 ) , R2 ( 3 ) , R3 ( 3 ) 

230  FOR  I1=-MAX  TO  MAX 

240  FOR  J1=-MAX  TO  MAX 

250  FOR  K1=-MAX  TO  MAX 

270  FOR  L=1  TO  3 

280  R1  (L)=I1*B(1,L)+J1*B(2,L)+K1*B(3,L) 

290  NEXT  L 

300  D1=R1 ( 1 ) *R1 ( 1 ) +R1 ( 2 ) *R1 ( 2 ) +R1 ( 3 ) *R1 ( 3 ) 

310  IF  (D1=0)  OR  (D1>M2)  THEN  710 

320  FOR  I 2= -MAX  TO  MAX 

330  FOR  J2=-MAX  TO  MAX 

340  FOR  K2=-MAX  TO  MAX 

350  FOR  L=1  TO  3 

360  R2(L)=I2*B(1,L)+J2*B(2,L)+K2*B(3,L) 

370  NEXT  L 

380  D2=R2  ( 1 )  *R2  ( 1 )  +R2  (  2 )  *R2  ( 2 )  +R2  (  3  )  *R2  (  3  ) 

390  IF  (D2=0)  OR  (D2>M2)  THEN  700 

400  FOR  I3=-MAX  TO  MAX 

410  FOR  J3=-MAX  TO  MAX 

420  FOR  K3  =  -MAX  TO  MAX 

430  FOR  L=1  TO  3 

440  R3(L)=I3*B(1,L)+J3*B(2, L)+K3*B( 3 , L) 

FIGURE  31.  Program  Ising4f  for  Calculating 

Coordination  Number  for  a  4-cluster 
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450  NEXT  L 

460  D3=R3 ( 1 ) *R3 ( 1 ) +R3 ( 2 ) *R3 ( 2 ) +R3 ( 3 ) *R3 ( 3 ) 

470  IF  (D3=0)  OR  (D3>M2)  THEN  690 
480  D12  =  0 
490  FOR  L=1  TO  3 

500  R12=R1(L)-R2(L) : D12=D12+R12*R12 
510  NEXT  L 

520  IF  (D12=0)  OR  (D12>M2)  THEN  690 
530  D13=0 
540  FOR  L=1  TO  3 

550  R13=R1(L)-R3(L) : D13=D13+R13*R13 
560  NEXT  L 

570  IF  (D13=0)  OR  (D13>M2)  THEN  690 
580  D23=0 
590  FOR  L=1  TO  3 

600  R23=R2 ( L ) -R3 ( L ) : D23=D23+R23*R23 
610  NEXT  L 

620  IF  (D23=0 )  OR  (D23>M2)  THEN  690 
630  REM  BUBBLE  SORT 
635  GOSUB  6000 
690  NEXT  K3, J3, 13 
700  NEXT  K2  , J2 , 12 
710  NEXT  Kl.Jl.Il 
720  PRINT  #1, "Dl,  D2 , 

730  FOR  Dl=l  TO  M2 
740  FOR  D2=l  TO  M2 
750  FOR  D3=l  TO  M2 
752  FOR  D12=l  TO  M2 
754  FOR  D13=l  TO  M2 
756  FOR  D23=l  TO  M2 
760  PRINT  #1,D1;"  M;D2;M  M;D3;"  ";D12;"  ";D13;" 

D23 ; ”  0(01,02,03,012,013,023) 

765  PRINT  Dl; "  ";D2;"  "jDS;"  ";D12;"  ";D13;M  ” ; 

D23 ; "  " ; C(D1 , D2 , D3 , D12 , D13 , D23 ) 

770  NEXT  023,013,012,03,02,01 
780  CLOSE  #1 

6000  D(1)=D1 :D(2)=D2:D(3)=D3:D(4)=D12:D(5)=D13: 

D( 6 ) =D23 : P1=D1 

6010  M=6 

6020  FOR  G=M-1  TO  1  STEP  -1 
6030  FOR  F=1  TO  G 

6040  IF  D(F)>D(F+1)  THEN  SWAP  D(F),D(F+1) 

6050  NEXT  F 
6060  NEXT  G 

6070  D1=D(1) :D2=D(2) :D3=D(3) :D12=D(4) :D13=D(5) :D23=D(6) 
6080  C (Dl , D2 , D3 , D12 , Dl 3 , D23 ) =C ( Dl , D2 , D3 , D12 , Dl 3 , D23 ) +1 
6085  D1=P1 
6090  RETURN 


D3 ,  D12 ,  D13 ,  D23 

Z(D1,D2,D3,D12,D13,D23)" 


FIGURE  31.  continued 
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APPENDIX  E 


PROGRAM  LISTINGS 

There  are  three  programs  listed  in  this  appendix. 

Fig.  32  is  a  listing  of  SIM4  used  to  obtain  the  equilibrium 
(minimum)  free  energy  for  the  fourth  approximation  which  is 
the  sum  of  the  first,  second,  third  and  fourth  term  of  the 
Morita  expansion.  The  minimization  program  for  lower  order 
approximations  can  be  obtained  by  deleting  the  equations  for 
the  higher  order  clusters.  The  Simplex  minimization  that  is 
used  here  is  described  in  reference  6. 

Fig.  33  is  a  listing  of  HIDDEN6  that  is  used  for 
making  the  contour  plots  shown  in  the  text.  This  program  is 
described  in  reference  [24-26] . 

The  program  GRID6  listed  in  Fig.  34  was  used  to 
calculate  the  entropy  values  in  the  third  and  fourth 
approximations  that  were  plotted  as  a  function  of  x,y.  . 
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10  REM  ** PROGRAM  NAME  IS  ’ SIM4 . BAS’ ,  BECAUSE  IT  DOES  A 
’SIMPLEX’  MINIMIZATION 
20  REM  **ON  THE  4-CLUSTER  ISING  MODEL 

30  REM  **COPIED  FROM  ’ SIM9 . BAS ’  RUNNING  ON  HEN’S  PDP-11 
40  REM  **21-MAR-1985 
50  RANDOMIZE 

60  Z1=12:Z111=48:Z111111=48 
70  TR=2 
80  BR=0 
90  GOTO  330 

99  REM  **PUT  THE  FUNCTION  TO-BE-MINIMIZED  HERE... 

100  ON  ERROR  GOTO  299 
110  X=A(0) 

120  Y=A(1) 

122  Z=A(2) 

124  W=A( 3 ) 

130  Pl=(l+X)/2  :  P2=(l-X)/2 
140  R1=P1*L0G(P1)+P2*L0G(P2) 

150  P21=  ( 1+2*X+Y) /4 : P22= ( 1 -Y) /4 : P23= ( 1-2*X+Y) /4 
160  R2=P21*LOG(P21 )+2*P22*LOG(P22 )+P23*LOG(P23 ) 

162  P31=  ( 1+3*X+3*Y+Z ) /8 : P34= ( 1-3*X+3*Y-Z ) /8 

163  P32  =  ( 1+X-Y-Z ) /8 : P33= ( 1-X-Y+Z ) /8 

164  R3=P31*LOG(P31 ) +P34*LOG(P34 ) +3*P32*LOG(P32 ) 
+3*P33*LOG(P33 ) 

170  P41=(1+4*X+6*Y+4*Z+W)/16:P43=( 1+2*X-2*Z-W) /16 
172  P42=(1-4*X+6*Y-4*Z+W)/16:P44=( 1-2*X+2*Z-W)/16 
174  P45=(1-2*Y+W)/16 

176  R4=P41*LOG(P41 )+P42*LOG(P42)+4*(P43*LOG(P43) 
+P44*LOG(P44) )+6*P45*LOG(P45) 

180  R=Rl+( Zl/2 ) *(R2-2*R1 ) +( Zlll/6 ) *(R3-3*R2+3*R1 ) 

182  R=R+(Z111111/24)*(R4-4*R3+6*R2-4*R1) 

190  F=-BR*X- ( Zl/2 ) *Y+TR*R 

298  ON  ERROR  GOTO  0 : RETURN 

299  F=lE+35: RESUME  298 
330  N=4 

340  11=0 

350  DIM  X(5,5) ,M(5) ,V(5) ,R(5) ,E(5) ,C(5) ,A(5) 

360  FOR  1=0  TO  N-l 
370  READ  M( I ) , V( I ) 

380  NEXT  I 

390  FOR  J=0  TO  N 

400  FOR  1=0  TO  N-l 

410  X(I,J)=M(I)+RND*V(I) 

420  NEXT  I 
430  NEXT  J 

440  F7=lE+37  :  F9=-F7 
450  V7=-l  :  V9=-l 
460  FOR  J=0  TO  N 


FIGURE  32.  Program  SIM4 
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470  FOR  1=0  TO  N-l 
480  A( I ) =X( I , J ) 

490  NEXT  I 

500  GOSUB  100 : X( N , J ) =F 
510  IF  F <F7  THEN  F7=F  :  V7=J 
520  IF  F>F9  THEN  F9=F  :  V9=J 
530  NEXT  J 

540  IF  V7  =  -l  THEN  PRINT"F-MIN  NOT  FOUND"  : STOP 

550  IF  V9=-l  THEN  PRINT "F-MAX  NOT  FOUND”  : STOP 

560  IF  V7=V9  THEN  PRINT"F-MIN  =  F-MAX !!!": STOP 

570  PRINT" ITTERATION  II : 11=11+1 
580  FOR  J=0  TO  N 
590  PRINT  J;"  "; 

600  FOR  1=0  TO  N 
610  PRINT  X( I , J) ; "  "  ; 

620  NEXT  I 
630  PRINT 
640  NEXT  J 
650  PRINT 

660  FOR  1=0  TO  N-l 
670  S=0 

680  FOR  J=0  TO  N 

690  IF  J<>V9  THEN  S=S+X(I,J) 

700  NEXT  J 

710  M(I)=S/N 

720  V(I)=X(I,V9)-M(I) 

730  R( I ) =M( I ) -V( I ) 

740  A( I ) =R( I ) 

750  NEXT  I 

760  GOSUB  100  :F1=F 

770  IF  F1<=F7  THEN  910 

780  IF  F1<=F9  THEN  1020 

790  FOR  1=0  TO  N-l 

800  C(I)=M(I)+V(I)/2  :  A(I)=C(I) 

810  NEXT  I 

820  GOSUB  100  :  F2=F 
830  IF  F2<=F9  THEN  1070 
840  FOR  1=0  TO  N-l 
850  FOR  J=0  TO  N 

860  IF  JOV7  THEN  X(  I ,  J)  =  (X(  I ,  J)+X(  I ,  V7  )  )/2 
870  NEXT  J 
880  NEXT  I 

890  PRINT "CONTRACT  ENTIRE  SIMPLEX” 

900  GOTO  1120 

910  FOR  1=0  TO  N-l 

920  E(I)=M(I)-2*V(I) : A( I ) =E ( I ) 

930  NEXT  I 
940  GOSUB  100:F3=F 
950  IF  F3<=F7  THEN  970 
960  GOTO  1020 


FIGURE  32  continued 
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970  FOR  1=0  TO  N-l 
980  X( I , V9 ) =E( I ) 

990  NEXT  I 

1000  PRINT " ACCEPT  EXPANDED  VERTEX" 
1010  GOTO  1120 
1020  FOR  1=0  TO  N-l 
1030  X( I , V9 ) =R( I ) 

1040  NEXT  I 

1050  PRINT "ACCEPT  REFLECTED  VERTEX" 
1060  GOTO  1120 
1070  FOR  1=0  TO  N-l 
1080  X( I , V9 ) =C ( I ) 

1090  NEXT  I 

1100  PRINT "ACCEPT  CONTRACTED  VERTEX" 

1110  GOTO  1120 

1120  GOTO  440 

1130  DATA  .0, .3 

1140  DATA  .0, .3 

1150  DATA  .0, .3 

1160  DATA  .0, .3 


FIGURE  32  continued 
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1  ’From  NIBBLE/Vol  4/No.  8/pp  61-71 

2  ’Xfered  "to  Jack’s  COMPAQ  on  3-27-86  by  HL 

3  ’ 

4  CMAX=639:RMAX=199 

5  Zl=12 

6  TR=10 

7  BR-0 

90  DIM  H(CMAX) ,L(CMAX) 

100  MM=1E+10 : UH=-MM: UL=MM: VH=-MM: VL=MM 
110  FOR  1=0  TO  CMAX : L ( I ) =RMAX : NEXT 
120  XE=2/2. 5:YE=l/2. 5: ZE=l/2. 5 

130  S1=XE*XE+YE*YE:S2=SQR(S1) : S3=SQR( S1+ZE*ZE ) : S4=l/ ( S2*S3 ) 
140  M=20 : N=20 

150  DIM  X(M) ,Y(N) ,R(M,N, 1) 

160  XL=-1:XH=1:YL=-1:YH=1 
170  DX= (XH-XL ) /M: DY= ( YH-YL) /N 
180  X0=XH : IF  XE<0  THEN  DX=-DX:X0=XL 
190  Y0=YH: IF  YE<0  THEN  DY=-DY : Y0=YL 
200  CX= ( CMAX+1 )/2 : CY= (RMAX+1 ) /2 
210  FOR  1=0  TO  M : X ( I ) =X0 - I *DX : NEXT 
220  FOR  J=0  TO  N : Y ( J ) = Y0 -  J*DY : NEXT 
230  FOR  1=0  TO  M: FOR  J=0  TO  N 


240  X=X(I):Y=Y(J) 

250  ’ 

260  ’  Replacable  function 
270  ’  Z=F(X, Y) 

272  ’ 

278  Pl=(l+X)/2:P2=(l-X)/2 

280  IF  P1<=0  OR  P2<=0  THEN  PHI=0:GOTO  300 


283  SSl=Pl*LOG(Pl ) +P2*LOG(P2 ) 

284  P21= ( l+2*X+Y)/4 : P22= ( 1-Y) /4 : P23= ( 1-2*X+Y) / 4 

285  IF  P21<=0  OR  P22<=0  OR  P23<=0  THEN  PHI=0:GOTO  300 

286  SS2=P21*LOG(P21 ) +2*P22*LOG(P22 )+P23*LOG(P23 ) 

288  PHI=-BR*X-(Z1/2)*Y+TR*(SS1+(Z1/2)*(SS2-2*SS1 ) ) 


300  IF  PHI>2  THEN  PHI=2 
310  Z=PHI/20+ . 5 
320  GOSOB  890 

330  R ( I , J , 0 ) =U : R( I , J , 1 ) = V : GOSUB  1420 
340  NEXT  J:  PRINT  " I  =  "  ;  I ,  ,,MAX="  ;  M:  NEXT  I 
350  ’S  is  the  scale  factor 
360  S=MM: IF  UL=0  THEN  380 

370  S0=CMAX/( 2 . 1*ABS(UL) ) : IF  S0<S  THEN  S=S0 
380  IF  UH=0  THEN  400 

390  S0=CMAX/(2. 1*ABS(UH) ) : IF  S0<S  THEN  S=S0 
400  IF  VL=0  THEN  420 

410  S0=RMAX/(2. 1*ABS(VL) ) : IF  S0<S  THEN  S=S0 
420  IF  VH=0  THEN  450 

430  S0=RMAX/( 2 . 1*ABS( VH) ) : IF  S0<S  THEN  S=S0 
440  ’Locate  in  HGR2  coordinates 
450  FOR  1=0  TO  M: FOR  J=0  TO  N 


FIGURE  33.  PROGRAM  HIDDEN6 
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460  R(I, J,0)=INT(S*R(I, J,0)+CX) :R(I, J, 1 )=INT(S*R( I , J, 1)+CY) 

470  NEXT  J : NEXT  I 

480  ’Start  graphics 

490  KEY  OFF: SCREEN  2 

492  CLS 

500  IF  ABS ( XE ) < ABS ( YE )  THEN  670 
510  FOR  1=0  TO  M 
520  ’Draw  fixed  X-lines 
530  FOR  J=1  TO  N 

540  U1=R(I, J-1,0) :V1=R(I, J-l, 1) :U2=R(I, J,0) :V2=R(I, J, 1) 

550  GOSUB  950: ’Test  visibility,  and  plot... 

560  GOSUB  1220: ’Update  H,L  arraws 

570  NEXT  J 

580  IF  I=M  THEN  650 

590  ’Draw  fixed  Y-line  segments 

600  FOR  J=0  TO  N 

610  U1=R(I, J,0) : V1=R(I, J, 1) : U2=R( 1+1 , J, 0 ) : V2=R( 1+1 , J, 1 ) 

620  GOSUB  950: ’Test  visibility,  and  plot... 

630  GOSUB  1220: ’Update  H, L  arrays 

640  NEXT  J 

650  NEXT  I 

660  GOTO  820 

670  FOR  J=0  TO  N 

680  ’Draw  fixed  Y-lines 

690  FOR  1=1  TO  M 

700  U1=R(I-1, J,0) :V1=R(I-1, J, 1) :U2=R(I, J,0) : V2=R(I, J, 1) 

710  GOSUB  950 

720  GOSUB  1220 

730  NEXT  I 

740  IF  J=N  THEN  810 

750  ’Draw  fixed  X-line  segments 

760  FOR  1=0  TO  M 

770  U1=R( I , J , 0 ) :V1=R(I, J, 1) : U2=R( I , J+l , 0 ) : V2=R( I , J+l , 1 ) 

780  GOSUB  950 
790  GOSUB  1220 
800  NEXT  I 
810  NEXT  J 

820  PRINT  CHR$ ( 7 ) : ’ Here  is  a  chance  to  get  printed  output. 
830  ’ 

840  ’ 

850  END 
860  ’ 

870  ’Transformation  subroutine 
880  ’ 

890  U= (XE*Y-YE*X) /S2 

900  V= ( ZE* (X*XE+Y*YE ) -S1*Z ) *S4 

910  RETURN 

920  ’ 

930  ’Wright’s  algorithm!!! 

940  ’ 

950  T1=0:T2=0:G1=0:G2=0 


FIGURE  33.  continued 
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1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 
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1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 
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1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 


LINE(U1,V1)-(U2JV2) : RETURN 
LINE(U1 , VI ) - (U2 , V2 ) : RETURN 
RETURN 


960  IF  V1>=H(U1 )  THEN  Tl=l 
970  IF  V2>=H(U2 )  THEN  T2=l 
980  IF  V1<=L(U1 )  THEN  Gl=l 
990  IF  V2<=L(U2)  THEN  G2=l 
1000  IF  Tl=l  AND  T2=l  THEN 
IF  Gl=l  AND  G2=l  THEN 
IF  T1+T2+G1+G2=0  THEN 
GOSUB  1370 
IF  KM=KX  THEN  1160 
F1=0 : F2=0 
FOR  K=KM  TO  KX 
VK=VM+(VX-VM)*(K-KM)/(KX-KM) 

IF  VK>H(K)  OR  VK<L(K)  THEN  U1=K: V1=VK: Fl=l : K=KX 
NEXT 

FOR  K=KX  TO  KM  STEP  -1 
VK= VM+ ( VX- VM ) * ( K-KM ) / ( KX-KM ) 

IF  VK>H(K)  OR  VK<L( K)  THEN  U2=K: V2=VK: F2=l : K=KM 
NEXT 

THEN  LINE(U1 ,V1)-(U2,V2) 


IF  Fl=l  AND  F2=l 
RETURN 

IF  VX>H(U1 )  THEN 
IF  VM<L(U1 )  THEN 
RETURN 


LINE(U1 ,H(U1) ) - (U1 , VX) : RETURN 
LINE (U1 ,L(U1))-(U1,VM) 


’Update  H,L  arrays 


H(U1)=V1 
H(U2 ) =V2 
L(U1)=V1 
L(U2 ) =V2 
THEN  RETURN 


IF  V1>H(U1 )  THEN 
IF  V2>H(U2 )  THEN 
IF  V1<L(U1)  THEN 
IF  V2<L(U2)  THEN 
IF  ABS(U1-U2 ) <1 
GOSUB  1370 
FOR  K=KM+1  TO  KX-1 
VK= VM+ ( VX- VM ) * ( K-KM ) / ( KX-KM ) 

IF  VK>H(K)  THEN  H(K)=VK 
IF  VK<L (K)  THEN  L(K)=VK 
NEXT  K 

RETURN 

> 

’Find  left-most  point  of  the  line 


KM=U1 : KX=U2 : VM=V1 : VX=V 2 : IF 
VM=V2: VX=V1: RETURN 
RETURN 


KM>KX  THEN  KM=U2:KX=U1: 


’Find  extreme  values  in  U,V  coordinates  before  scaling 


IF 

IF 

IF 

IF 


U>UH 
U<UL 
V>  VH 
V<VL 


THEN 

THEN 

THEN 

THEN 


UH=U 

UL=U 

VH=V 

VL=V 


RETURN 


FIGURE  33  continued 
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10  Zl=12 
20  Zlll=48 
30  Zllllll=48 

100  OPEN  "lptl:"  FOR  OUTPUT  AS  # 1 
102  LG2=LOG( 2 ) 

105  FOR  W= . 8  TO  0  STEP  -.2 

106  FOR  Z=.8  TO  0  STEP  -.2 
108  PRINT  #1 , "Z=” ; Z; "w=" ;W; 

110  PRINT  #1,"  Y-AXIS" 

112  PRINT  #1,"  M; 

120  FOR  Y=-10  TO  10 

122  PRINT  #1,  USING  Y; 

124  NEXT  Y 
126  PRINT  #1, 

200  FOR  X=1  TO  -1  STEP  -.1 

202  PRINT  #1 , USING  X; : PRINT  #1,” 

210  FOR  Y=-l  TO  1.05  STEP  .1 
220  ON  ERROR  GOTO  800 
300  Pl=(l+X)/2 
310  P2=(l-X)/2 

320  R1=P1*L0G(P1)+P2*L0G(P2) 

330  P21=(l+2*X+Y)/4 
340  P22=(l-Y)/4 
350  P23=(l-2*X+Y)/4 

360  R2=P21*LOG(P21 ) +2*P22*LOG(P22 ) +P23*LOG(P23 ) 

400  P31=(l+3*X+3*Y+Z)/8 
410  P32= ( 1+X-Y-Z ) /8 
420  P33= ( 1-X-Y+Z ) /8 
430  P34=(l-3*X+3*Y-Z)/8 

440  R3=P31*LOG(P31)+3*P32*LOG(P32)+3*P33*LOG(P33) 
+P34*LOG(P34 ) 

445  P41= ( 1+4*X+6*Y+4*Z+W)/16 : P43= ( 1+2*X-2*Z-W) /16 
450  P42=(1-4*X+6*Y-4*Z+W)/16:P44=(1-2*X+2*Z-W)/16 
453  P45=(1-2*Y+W)/16 

455  R44=P41*LOG(P41 )+P42*LOG(P42 ) +4* (P43*LOG(P43 ) 
+P44*LOG(P44) )+6*P45*LOG(P45) 

460  R4=R1+(Z1/2)*(R2-2*R1)+(Z111/6)*(R3-3*R2+3*R1) 
465  R4=R4+( Zll 1111/24 )*(R44-4*R3+6*R2-4*R1 ) 

470  R4=-R4/(LG2) 

600  R=INT ( 100*R4+ . 5 ) 

610  PRINT  #1 , USING  "###" ;R; 

700  NEXT  Y 
710  PRINT  #1, 

720  NEXT  X 
730  NEXT  Z 
740  NEXT  W 

799  END 

800  PRINT  #1,”  **"; 

810  RESUME  700 


FIGURE  34.  PROGRAM  GRID6 
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APPENDIX  F 


UNITS  and  NOTATION 

The  following  "reduced"  units  are  used  in  the  text. 

5  =  "reduced"  free  energy  F  =  (F/NJ) 

Br  =  "reduced"  magnetic  field  B  =  (J-'oB/J) 

Tr  =  "reduced"  temperature  T  =  (kBT/J) 

J  =  the  value  of  the  exchange  integral  =  (kBTc/Z). 

For  a  Curie  temperature  =  Tc  =  1056  °K  (Fe)  and  Z  =  12  for  an 
FCC  lattice:  J  =  7.56  (10-3)e.v. 

S(  Joules/°K) 

Entropy:  - -  S  in  "bits" 

In  2 

In  the  equations  two  notations  have  been  used  to 
represent  the  same  quantites. 

either  x 

<P>  = 

or  x<  i  > 

either  yi  j 

<PP>  = 

or  x<  2  >  i  j 

either  zi  j  k 

<PPP>  = 

or  x<  3  >  i  j  k 
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APPENDIX  G 


TABLE  2 

EQUILIBRIUM  FREE  ENERGY  <5 
Br  -  0 


Tr 

Approximation 

First  J  Second 

|  Third(^) 

|  Fourth 

2 

-6.000012 

- 6 . xxxxxx 

-5.991287 

-5.734586 

4 

-6.010053 

-6.010166 

-6.043771 

.  -6.009864 

6 

-6.118027 

-6.122647 

-6.118978 

i  -6.123150 

8 

-6.466731 

-6.500590 

-6.503969 

-6.449281 

10 

-7.172468 

-7.298410 

-7.319456 

-7.375871 

12 

-8.317766 

-8.567490 

-8.599178 

-8.653156 

14 

-9.704059 

:  -9.918160 

xxxxxxxxx 

-9.975911 

16 

-11.09035 

-11.27774 

-11.21078 

-11.31896 

18 

-12.47665 

-12.64323 

-12.65670 

-12.67498 

20 

-13.86294 

-14.01289 

-14.02367 

-14.03782 

3 

- 

-6.001010 

- 

— 

5 

- 

-6.043730 

- 

- 

7 

- 

-6.268225 

-6.269251 

- 

9 

- 

-6.838270 

-6.847406 

- 

9.5 

- 

-7.052020 

- 

- 

10.25 

- 

-7.434480 

- 

- 

10.5 

- 

-7.579450 

- 

- 

10.75 

- 

-7.733550 

- 

- 

10.85 

- 

-7.797817 

- 

- 

10.9 

- 

-7.805160 

- 

- 

10.95 

- 

-7.863579 

- 

- 

11 

- 

-7.896960 

-7.935181 

- 

13 

- 

-9.241450 

-9.268160 

- 

15 

- 

-10.59706 

-10.61683 

- 

17 

- 

-11.95988 

-11.97504 

- 

* 
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TABLE  3. 


SOME  CORRELATION  COEFFICIENTS 
Br  -  0 


Cluster 

Correlation  Coefficients 

X 

X 

1  y  1  z  1 

w 

Tr  =  4 


1- c 

2- c 

3- c(*) 

4- c 

.9948288 
.9947500 
. 9651874 
.9945275  } 

.9895800 
. 9330112 
.9891984 

. 9030889 
. 9839292 

.  9787838 

03 

It 

u* 

tr* 

1- c 

2- c 

3- c(S) 

4- c 

. 8581825 
. 8325000 
. 8276450 
. 7032843 

.7070200 

.7017475 

.5307391 

. 6039199 
.4117586 

.  3277481 

J  Tr  =  10 

1- c 

2- c 

3- o($) 

4- c 

.6588523- 
.5350000 
.4924544 
.  1001604 

.3400000 
. 3144164 
.2064547 

. 2087749 
.0430753 

■ 

.08437848 

Tr  =  16 


1-c 

1 . 089163E-3 

- 

- 

- 

2-c 

-1 . 298148E-3 

6 . 244406E-2 

- 

- 

3-c(%) 

-.1772807 

9. 670085E-2 

-2. 161592E-2 

- 

4-c 

9 . 9860E-5 

8 . 6007E-2 

-1. 988E-4 

1 . 8118E-2 

t 
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